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Abstract 

The  thesis  develops  and  demonstrates  methods  of  classifying  ocean  processes  using 
an  underwater  moving  platform  such  as  an  Autonomous  Underwater  Vehicle  (AUV). 
The  “mingled  spectrum  principle”  is  established  which  concisely  relates  observations 
from  a  moving  platform  to  the  frequency-wavenumber  spectrum  of  the  ocean  pro¬ 
cess.  It  clearly  reveals  the  role  of  the  AUV  speed  in  mingling  temporal  and  spatial 
information.  For  classifying  different  processes,  an  AUV  is  not  only  able  to  jointly 
utilize  the  time-space  information,  but  also  at  a  tunable  proportion  by  adjusting 
its  cruise  speed.  In  this  respect,  AUVs  are  advantageous  compared  with  traditional 
oceanographic  platforms. 

Based  on  the  mingled  spectrum  principle,  a  parametric  tool  for  designing  an  AUV- 
based  spectral  classifier  is  developed.  An  AUV’s  controllable  speed  tunes  the  separa¬ 
bility  between  the  mingled  spectra  of  different  processes.  This  property  is  the  key  to 
optimizing  the  classifier’s  performance. 

As  a  case  study,  AUV-based  classification  is  applied  to  distinguish  ocean  convec¬ 
tion  from  internal  waves.  The  mingled  spectrum  templates  are  derived  from  the  MIT 
Ocean  Convection  Model  and  the  Garrett-Munk  internal  wave  spectrum  model.  To 
allow  for  mismatch  between  modeled  templates  and  real  measurements,  the  AUV- 
based  classifier  is  designed  to  be  robust  to  parameter  uncertainties.  By  simulation 


tests  on  the  classifier,  it  is  demonstrated  that  at  a  higher  AUV  speed,  convection’s 
distinct  spatial  feature  is  highlighted  to  the  advantage  of  classification. 

Experimental  data  are  used  to  test  the  AUV-based  classifier.  An  AUV-borne  flow 
measurement  system  is  designed  and  built,  using  an  Acoustic  Doppler  Velocimeter 
(ADV).  The  system  is  calibrated  in  a  high-precision  tow  tank.  In  February  1998,  the 
AUV  acquired  field  data  of  flow  velocity  in  the  Labrador  Sea  Convection  Experiment. 
The  Earth-referenced  vertical  flow  velocity  is  extracted  from  the  raw  measurements. 
The  classification  test  result  detects  convection’s  occurrence,  a  finding  supported  by 
more  traditional  oceanographic  analyses  and  observations.  The  thesis  work  provides 
an  important  foundation  for  future  work  in  autonomous  detection  and  sampling  of 
oceanographic  processes. 
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Chapter  1 


Introduction 

1.1  Motivation 

One  of  the  most  challenging  tasks  in  observing  and  studying  the  ocean’s  temporal  and 
spatial  variability  is  to  identify  the  underlying  ocean  processes.  The  thesis  develops 
and  demonstrates  methods  of  classifying  ocean  processes  using  observations  from  an 
Autonomous  Underwater  Vehicle  (AUV).  Automated  classification  will  also  enable 
adaptive  sampling  [1]  and  other  refined- monitoring  measures  [5],  [6]. 

Eulerian  and  Lagrangian  platforms  are  representative  of  traditional  oceanographic 
monitoring  tools  [7].  An  Eulerian  platform  is  fixed  in  location,  providing  time  series 
records  of  measured  quantities.  Moored  current  meters  and  Conductivity-Temperature- 
Depth  (CTD)  sensors  have  become  a  routine  in  oceanographic  monitoring.  A  La¬ 
grangian  platform,  on  the  other  hand,  drifts  with  the  current  flow.  The  path  followed 
by  Lagrangian  drifters  reveals  the  current  flow  history  [8],  [9].  By  tracking  Lagrangian 
platforms  acoustically  (e.g.  SOFAR  drifters)  or  by  satellite  (e.g.,  the  ARGOS  system) 
for  surface  floats,  we  can  obtain  a  first-order  description  of  the  global  ocean  circula¬ 
tion.  Surface  drifters  and  many  subsurface  floats  are  considered  “quasi-Lagrangian” 
since  they  move  on  a  two-dimensional  plane  [7].  Improved  Lagrangian  floats  [10]  fol¬ 
low  the  three  dimensional  motion  of  water  parcels  through  density  matching  and  a 
high  drag  (we  will  see  example  results  in  Section  9.3). 

Classification  imposes  a  higher  level  of  requirement  than  monitoring.  Both  tem- 
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poral  and  spatial  features  should  be  utilized  to  optimize  classification.  Eulerian  and 
Lagrangian  platforms  have  inherent  limitations  in  this  respect.  Eulerian  measurement 
is  confined  to  a  fixed  location.  Although  a  mooring  may  sense  some  information  of 
the  field’s  spatial  variation  via  a  horizontal  advective  current,  this  kind  of  sensing  is 
uncontrolled  and  tends  to  be  ambiguous.  Deploying  an  array  of  moorings  can  add  in 
spatial  coverage,  but  high  cost  would  often  deter  dense  spatial  sampling  (please  see 
the  example  under  Item  1  in  Section  2.1).  A  Lagrangian  platform  drifts  with  zero 
relative  velocity  against  the  ambient  flow.  It  does  move,  but  its  motion  is  no  different 
from  the  advecting  current.  As  a  drifter  is  bound  to  a  tagged  parcel  of  water,  it 
has  hardly  any  chance  to  catch  sight  of  the  real  spatial  variation  of  the  field.  Since 
Eulerian  or  Lagrangian  platforms  have  limitations  in  providing  temporal  plus  spatial 
features  of  ocean  processes,  we  resort  to  moving  platforms  to  overcome  this  deficiency. 

A  towed  platform  is  tied  to  a  surface  ship.  This  type  of  platform  is  typically 
confined  to  a  depth  of  no  more  than  a  few  hundred  meters  [11].  A  larger  depth 
slows  the  tow  speed,  limits  maneuverability,  and  increases  the  cost  of  the  cable  and 
winch  system.  An  Autonomous  Underwater  Vehicle  (AUV)  [1]  is  an  unmanned, 
untethered  moving  platform.  The  Odyssey  IIB  AUVs  can  dive  to  the  full  ocean 
depth  in  most  places  [12].  An  AUV  is  neither  Eulerian  nor  Lagrangian,  but  cruises 
through  the  ocean  at  a  controllable  and  flexible  speed,  collecting  information  of  both 
time  and  space.  AUV  measurements  mingle  the  temporal  and  spatial  variations  of 
the  field.  Once  equipped  with  a  classification  capability,  an  AUV  has  the  promise 
of  autonomously  searching  for  oceanographic  processes  of  interest.  Hence  we  are 
encouraged  to  explore  how  to  utilize  its  controllable  mingling  of  temporal  and  spatial 
information  to  the  advantage  of  classification.  The  thesis  addresses  this  problem 
with  goals  of:  developing  a  parametric  tool  for  designing  an  AUV-based  classifier, 
and  demonstrating  the  AUV-based  classification  by  simulations  and  experiments. 


20 


1.2  Thesis  Outline 


Chapter  2  reviews  existing  research  in  relation  to  the  thesis  work.  An  AUV  is  com¬ 
pared  with  traditional  oceanographic  platforms:  Eulerian,  Lagrangian,  and  towed. 
Previous  formulations  of  the  Doppler  effect  on  spectrum  measurement  are  introduced. 
Another  area  reviewed  is  statistical  classification,  which  is  the  background  for  design¬ 
ing  the  AUV-based  classifier.  Theories  of  feature  extraction  and  linear  classifier  lead 
to  the  spectral  feature  classification  method  to  be  presented  in  Chapter  4. 

Chapter  3  develops  the  mingled  spectrum  principle.  Derivations  and  interpreta¬ 
tions  are  given.  The  comparison  with  previous  research  is  discussed.  The  presented 
formula  concisely  relates  observations  from  a  moving  platform  to  the  temporal-spatial 
spectrum  of  the  process  under  survey.  It  clearly  reveals  the  role  of  the  AUV  speed  in 
mingling  temporal  and  spatial  information.  The  mingled  spectrum  principle  lays  the 
theoretical  basis  for  AUV-based  classification. 

Chapter  4  presents  the  design  of  an  AUV-based  spectral  feature  classifier.  We  ap¬ 
ply  the  general  method  of  feature  extraction  to  the  classification  of  Power  Spectrum 
Density  (PSD)  estimates.  The  class  separability  metric  and  the  resultant  transfor¬ 
mation  vector  are  built  upon  the  statistics  of  PSD  estimates.  It  is  the  AUV  speed 
that  tunes  the  separability  of  the  mingled  spectra  of  different  processes.  This  fact  is 
the  key  to  optimizing  the  classifier’s  performance. 

Chapter  5  introduces  two  oceanographic  processes  to  demonstrate  AUV-based 
classification.  They  are  ocean  convection  and  internal  waves.  The  vertical  flow  ve¬ 
locity  is  the  measured  quantity  used  for  classification.  We  run  the  MIT  Convection 
Model  using  experimental  parameters.  The  temporal-spatial  spectrum  of  convective 
vertical  velocity  is  obtained  from  the  model  output.  For  internal  waves,  we  apply  the 
well-known  Garrett-Munk  model.  The  temporal-spatial  spectrum  of  internal  wave 
vertical  velocity  is  derived  from  the  original  formulation. 

Chapter  6  merges  the  three  preceding  chapters.  Given  the  temporal-spatial  spec¬ 
tra  of  ocean  convection  and  internal  waves,  the  mingled  spectra  observed  by  the  AUV 
are  found.  The  AUV-based  classifier  is  applied  to  distinguish  these  two  oceanographic 
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processes.  Model  data  are  used  in  the  simulations.  The  simulation  results  demon¬ 
strate  that  at  a  higher  vehicle  speed  the  classification  performance  is  better  since  the 
distinction  between  convection  and  internal  waves  is  highlighted. 

Chapter  7  adds  robustness  to  the  classifier.  Because  of  parameter  uncertainties, 
the  classifier  needs  to  be  robust  to  model  mismatch.  Model  parameter  uncertainties 
are  categorized  into  “local  uncertainty”  and  “global  uncertainty”,  referring  to  small 
perturbation  and  major  mismatch,  respectively.  The  effects  of  these  uncertainties  are 
to  increase  the  total  variance  of  PSD  estimates,  and  accordingly  suppress  the  ampli¬ 
tude  of  the  feature  transformation  vector.  The  classifier’s  performance  (probability  of 
detection  versus  probability  of  false  alarm)  is  thus  lowered,  but  with  a  gain  of  robust¬ 
ness  to  model  mismatch.  Through  progressive  work  from  Chapter  3  to  Chapter  7,  we 
have  developed  a  parametric  tool  for  designing  an  AUV-based  spectral  classifier. 

Chapter  8  proceeds  to  the  experimental  work.  We  built  an  AUV-borne  flow  veloc¬ 
ity  measurement  system  using  an  Acoustic  Doppler  Velocimeter  (ADV).  To  verify  the 
vehicle  hull’s  influence  on  flow  measurement,  we  conducted  a  calibration  experiment 
in  the  David  Taylor  Model  Basin  —  a  high-precision  tow  tank.  The  data  processing 
method  to  extract  the  Earth-referenced  vertical  flow  velocity  from  the  AUV’s  raw 
measurements  is  presented.  The  error  analysis  is  formulated. 

Chapter  9  brings  the  AUV  to  the  open  ocean.  In  January/February  1998,  the 
AUV-borne  flow  velocity  measurement  system  acquired  data  from  the  Labrador  Sea. 
Meteorological  and  hydrographic  conditions  for  convection,  as  well  as  Lagrangian 
float  observations  of  convection,  are  described  in  relation  to  AUV  missions.  The 
Earth-referenced  vertical  flow  velocity  is  extracted  from  the  AUV’s  measurements. 
The  experimental  data  are  used  to  test  the  AUV-based  classifier.  The  test  result  for 
the  250-m  depth  measurement  during  AUV  Mission  B9804107  detects  convection’s 
occurrence. 

Chapter  10  summarizes  the  thesis  and  proposes  future  work. 
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Chapter  2 

Review  of  Existing  Work 


2.1  Autonomous  Underwater  Vehicle  in  Compari¬ 
son  with  Traditional  Oceanographic  Platforms 


Figure  2.1:  An  Odyssey  IIB  AUV  being  recovered  after  operations  [1]. 

An  Autonomous  Underwater  Vehicle  (AUV)  is  a  mobile  instrument  platform  [1]:  An 
Odyssey  IIB  AUV  [12]  (designed  and  built  by  the  MIT  Sea  Grant  AUV  Laboratory) 
is  shown  in  Figure  2.1.  Each  vehicle  has  a  length  of  2.2  m,  and  a  diameter  of  0.6  m:  at 
its  largest  vertical  cross-section.  The  vehicle  has  an  outer  fairing  for  hydrodynamic 
stability  and  drag  reduction,  and  an  inner  fairing  for  structural  integrity.  These  struc- 
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tures  are  free-flooded,  except  for  two  17-inch-diameter  glass  spheres  which  provide 
the  dry  volume.  The  vehicles  operate  autonomously  with  its  on-board  computer  and 
battery  system,  running  various  sensors.  With  minimal  logistical  support,  they  can 
be  deployed  at  remote  sites,  off  ships  of  opportunity,  or  in  rough  seas.  In  the  past 
several  years,  Odyssey  IIB  AUVs  have  fulfilled  about  400  missions  under  the  Arctic 
ice,  over  the  Pacific  Ocean  ridge,  in  the  Haro  Strait  tidal  current,  off  the  coast  of  New 
Zealand,  at  the  polar  Labrador  Sea,  and  along  the  narrow  Monterey  Canyon. 

A  comparison  between  an  AUV  and  traditional  oceanographic  platforms  is  made 
in  Table  2.1.  Note  that  the  comparison  of  temporal-spatial  capability  is  from  the 
specific  perspective  of  process  classification,  since  good  classification  requires  jointly 
utilizing  temporal  and  spatial  information.  To  help  explain  Table  2.1,  we  briefly 
review  traditional  oceanographic  platforms  as  follows. 


Table  2.1:  Comparison  of  an  AUV  with  traditional  oceanographic  platforms 


Records  field’s 
temporal  variation 

Records  field’s 
spatial  variation 

Deep  ocean 
capacity 

Platform 

motion 

Autonomous 

AUV 

V 

V 

V 

Controlled 

V 

Mooring 

(Eulerian) 

V 

V 

Disturbed  by 
ambient  current 

V 

Drifter 

(Lagrangian) 

/ 

S 

V 

Follows 
current  flow 

V 

Towed 

platform 

V 

V 

V* 

Tied  to 
surface  ship 

T-  As  a  drifter  is  bound  to  a  tagged  parcel  of  water,  it  can  hardly  sense 
the  real  spatial  variation  of  the  field. 


*:  Typically  no  more  than  a  few  hundred  meters  [11]  deep.  A  larger 

depth  slows  the  tow  speed,  limits  maneuverability,  and  increases  the  cost  of 
the  cable  and  winch  system. 


1.  Eulerian  platform. 

An  Eulerian  platform  is  fixed  in  location,  providing  time  series  records  of  mea¬ 
sured  quantities  (such  fixed-location  measurement  is  named  after  L.  Euler  who 
first  formulated  the  fluid  motion  equations  in  a  fixed  frame  of  reference  [7]). 
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Moorings  are  the  most  common  Eulerian  platforms.  Moored  Conductivity- 
Temperature-Depth  (CTD)  sensors  [13]  and  current  meters  [14]  have  become 
a  routine  in  oceanographic  monitoring  [7],  [15].  While  being  able  to  provide 
long-term  observations,  moored  measurement  is  confined  to  one  fixed  location. 
Although  a  mooring  may  sense  some  information  of  the  field’s  spatial  variation 
via  a  horizontal  advective  current,  this  kind  of  sensing  is  uncontrolled  and  tends 
to  be  ambiguous.  Deploying  an  array  of  moorings  can  add  in  spatial  coverage, 
but  high  cost  would  often  deter  dense  spatial  sampling.  During  the  Northern 
Mediterranean  deep  convection  experiment  in  winter  1991/92,  a  triangular  array 
of  current  meter  moorings  were  deployed.  The  2-km  spacing  between  moorings 
was  so  large  that  the  plume  measurements  made  at  the  three  stations  were 
decorrelated.  Hence  the  three  stations’  measurements  could  not  be  combined 
into  a  joint  plume  analysis  [16]. 

2.  Lagrangian  platform. 

A  Lagrangian  platform  drifts  with  the  current  flow,  thus  can  be  considered 
following  a  tagged  parcel  of  water  (such  path-following  measurement  is  named 
after  J.  Lagrange  who  is  noted  for  his  early  work  on  fluid  dynamics  [7]).  Surface 
floats  can  be  tracked  by  satellites  (e.g.,  the  ARGOS  system)  while  subsurface 
drifters  are  tracked  acoustically,  such  as  SOFAR  (SOund  Fixing  And  Ranging) 
drifters  which  are  neutrally  buoyant  at  the  sound  channel  depth.  Tracking 
of  a  large  number  of  drifters  provides  a  first-order  description  of  the  ocean 
current  [17].  It  seems  that  a  Lagrangian  platform  obtains  both  temporal  and 
spatial  information  by  moving  through  the  field,  but  it  cannot  sense  the  real 
spatial  variation  of  the  field  because  it  is  bound  to  a  tagged  parcel  of  water. 
Launching  a  large  number  of  drifters  may  help  construct  a  temporal-spatial 
picture  of  the  current,  but  this  picture  should  be  taken  only  in  a  statistical 
mean  sense  because  of  variability  of  particle  trajectories  due  to  the  ocean’s 
large  Reynolds  number  [18]. 

3.  Towed  platform. 
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A  towed  platform  is  hauled  by  a  surface  ship.  It  utilizes  the  ship  speed  to 
carry  out  underwater  surveys.  Towed  measurement  has  been  used  for  numerous 
oceanographic  studies  [19],  [20],  including  on  internal  waves  [2],  A  recent  ex¬ 
ample  of  towed  platform  is  Seasoar  [11]  which  makes  CTD  profiles  in  the  upper 
350-m  water  column.  It  is  equipped  with  controllable  wings  to  enable  ascent  or 
descent  to  desired  depths.  In  the  respect  of  motion,  a  towed  platform  resembles 
an  AUV,  but  it  has  the  following  constraints  that  limit  its  applicability: 

•  It  is  typically  confined  to  a  depth  of  no  more  than  a  few  hundred  me¬ 
ters  [11].  A  larger  depth  would  add  loads  on  the  cable  and  wrench  system. 

•  Its  motion  is  essentially  affected  by  that  of  the  surface  vessel.  Accurate 
motion  control  is  hard  to  achieve. 

•  Its  operation  completely  relies  on  a  ship.  Sustained  human  attendance  and 
expensive  ship  time  make  such  operations  costly. 


2.2  Doppler  Effect  on  Measurement’s  Spectrum 

Internal  wave’s  power  spectrum  is  established  in  [2]  based  on  moored,  towed,  and 
dropped  measurements  in  the  ocean.  The  paper  considers  the  towed  spectrum  from 
the  perspective  of  Doppler  effect  on  elementary  waves.  Its  derivation  is  summarized 
as  follows,  with  some  variable  notations  changed  to  those  used  in  the  thesis. 

For  a  sensor  moving  in  the  positive  x1  direction  (in  a  horizontal  plane)  at  a  con¬ 
stant  speed  u,  one  has  xx  —  ut  where  t  is  time.  The  observed  phase  of  an  elementary 
wave  train  e-77  is  then 

7  =  kxxx  +  k2x2  cot  —  k2x2  -  (w  -  kxu)t  (2.1) 

where  to  is  the  angular  frequency;  kx,  k2  are  the  two  orthogonal  horizontal  wavenum¬ 
bers.  Denote  the  radial  wavenumber  as  k,  and  then  k 2  =  kf  +  k%.  In  [2],  it  is  assumed 
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that  the  field  is  horizontally  isotropic  so  that  k  suffices  for  the  description  of  spatial 
variation. 

Denote  the  frequency-of-encounter  as  a,  and  then 


a  =  \u)  —  k\u\  (2-2) 

Power  contributions  to  a  come  from  three  branches  (allowing  that  u  is  always 
positive): 

1.  k\  <  0,  co  >  k\u  (=>  ui  <  a):  negative  waves  (traveling  opposite  to  the  sensor) 

2.  ki  >  0,  to  >  k\U  (=>•  u>  >  cr):  positive  waves,  overtaking  the  sensor 

3.  k\  >  0,  ui  <  k\U :  the  sensor  overtaking  positive  waves 

Accordingly,  the  towed  spectrum  of  internal  wave’s  vertical  displacement  is  com¬ 
posed  of  three  parts  as  follows.  Note  that  due  to  isotropy,  wavenumber  integrations 
are  over  the  radial  wavenumber  k. 

Ft(a)  =^Tr~3M~3u~l  J  Z2(u)dulJ  [k2  — 

1  u 

POO  _ 

+  /  [k2  —  ( - )2]-2  E(k,u)dk 

J  u  —  cr  U 

u 

+  [  [ k 2  -  (— +  -)2]~^E(k,u)dk\ 

J  u'+cr  U  J 

u 

where  Ui  is  the  Coriolis  frequency  and  n  is  the  buoyancy  frequency  at  the  studied 
depth  ( n  >  cu*).  These  two  bounds  confine  the  internal  wave’s  frequency  range.  M  is  a 
scaling  wavenumber  corresponding  to  the  e-folding  depth  of  the  buoyancy  frequency. 
Z(u)  is  a  real  wave  function  of  the  vertical  displacement,  and  Z2(u>)  is  the  mean- 
square  averaged  over  modes.  E(k,u)  is  the  energy  density  spectrum.  Coefficient  u~l 
is  due  to  the  scaling  between  a  and  ki  as  in  Equation  (2.2). 


(— — —)2}  2 E(k,uj)dk 
u 

(2.3) 
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Equation  (2.3)  thus  gives  a  method  of  calculating  the  sensor-recorded  spectrum 
from  an  isotropic  spectrum  of  the  original  field,  as  Z2(cu)E(k,  ui)  can  be  re¬ 

garded  as  the  spatial-temporal  spectrum  Sra(nai(k,uj).  By  removing  frequency  bound 
limitations  associated  with  specific  processes,  and  adopting  notations  Sraciiai{k,u )  for 
the  original  spatial-temporal  spectrum  and  S(a)  for  the  sensor- recorded  spectrum, 
Equation  (2.3)  can  be  generalized  to  (although  [2]  has  not  explicitly  made  this  ef¬ 
fort): 


S(a) 


j  poo  r  poo 

=~  d^\  [k2-{ - -)2}-^Sradial(h,uj)dk 

UJ- OO  u 

u 

POO 

+  [k2  -  ( - )2]~iSradiai(k,u))dk 

J  U 

U 

+  ^+tr^2  ~  ('  ^  )2]~^>S'radia;(^»a;)^| 

u  ' 


(2.4) 


Integral  I 
Integral  II 


Figure  2.2:  Integration  regions  in  calculation  of  the  Doppler-shifted  spectrum  (based 
on  Figure  8  of  [2]). 

In  Equation  (2.4),  the  three  integrals  enclosed  in  curly  braces  correspond  to  the 
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three  branches  listed  below  Equation  (2.2).  The  integrations  are  illustrated  in  Fig¬ 
ure  2.2,  where  the  three  integral  components  in  Equation  (2.4)  are  marked  by  region 
“I”,  “II”,  and  “III”,  respectively.  In  Figure  2.2,  integral  region  I  is  to  the  right  of 
line  (1);  region  II  to  the  right  of  line  (2);  and  region  III  to  the  right  of  line  (3).  Line 
u)  =  a  bisects  regions  I  and  II.  When  u  =  0  (moored  measurement),  lines  (1)  and  (2) 
coalesce  into  the  horizontal  line  im  =  a.  When  u  — »  oo,  lines  (1)  and  (3)  approach  the 
vertical  line  k  =  afu.  Note  that  region  III  overlaps  with  regions  I  and  II. 

In  [2],  Equation  (2.3)  is  called  “Doppler-shifted  spectrum”.  In  Subsection  3.1.3,  we 
will  discuss  the  differences  between  the  Doppler-shifted  spectrum  method  introduced 
herein  and  the  mingled  spectrum  principle  to  be  presented  in  Chapter  3. 

2.3  Statistical  Classification 

The  objective  of  classification  is  to  determine  which  class  a  given  sample  set  belongs 
to  [21],  [3].  Through  a  measurement  process,  we  obtain  an  observation  vector.  By 
a  decision  rule  we  assign  the  observation  vector  to  one  of  the  given  classes.  If  the 
conditional  probability  density  function  for  each  class  is  known,  the  classification 
problem  becomes  one  of  statistical  hypothesis  testing. 

2.3.1  Bayes  Decision  Rule 

In  this  thesis,  we  study  two-class  problems.  Denote  the  classes  as  H\  and  H2,  and 
the  observation  vector  as  X.  The  conditional  probability  density  functions  of  X  are 


Under  Hx  :  p(X\Hx ) 

(2.5) 

Under  H2  :  p(X\H2) 

(2.6) 

The  ratio  is  termed  the  “likelihood  ratio” . 

The  a  priori  probabilities  of  the  two  classes  are  denoted 
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PrlHi)  =  P, 

Pr{H2}  =  P2 


(2.7) 

(2.8) 


The  task  is  to  determine  whether  X  belongs  to  Hx  or  H2.  A  wrong  decision  is 
penalized  by  an  associated  cost.  The  classification  decision  should  be  made  such  that 
the  total  average  cost  (in  the  statistical  sense)  is  minimized.  The  cost  associated  with 
each  decision  is  expressed  as 


Cij  =  cost  of  deciding  X  e  Hi  when  X  e  Hj  (2.9) 

Usually  a  correct  decision  incurs  no  cost,  i.e.,  cn  =  c2 2  =  0. 

The  criterion  in  devising  the  Bayes  decision  rule  is  to  minimize  the  total  average 
cost,  also  called  the  “Bayes  cost”: 


R(f)  =  E[C(f(X),H)] 

=  cuPr{H x,  Hx }  +  c12Pr{#x,  H2}  +  c21Pr{H2,  Hx}  +  c22Pr{H2 ,  H2 ^ 

where  C(Hi,Hj )  =  (as  defined  in  Equation  (2.9))  is  the  cost  function;  /  denotes 
the  decision  rule;  f(X)  takes  on  either  Hx  or  H2  ('is  added  to  denote  the  decision, 
which  may  or  may  not  be  correct);  H  denotes  the  truth,  taking  on  either  Hx  or  H2. 
Note  that  the  expectation  is  executed  on  both  X  and  H. 

It  can  be  deduced  that  the  “Bayes  decision  rule”  [21],  [3]  minimizes  R(f).  The 
rule  is  expressed  by 


P{X\H2)  A2  (c2i  -  Ci x) Pi 
P{X\Hi)  g  (ci2  -  c22)P2 
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h2 

where  ^  means  that  when  “>”  holds,  we  decide  iL2;  otherwise  we  decide  H\.  The 
Hi 

above  rule  is  also  called  a  “Likelihood  Ratio  Test  (LRT)”. 


2.3.2  Bayes  Error  and  Its  Upper  Bound 

When  Cij  =  1  —  ( 5^  is  the  Kronecker  Delta  notation),  i.e.,  correct  decisions  carry 
no  cost  while  incorrect  decisions  induce  a  cost  of  one,  then  the  Bayes  cost  in  Equa¬ 
tion  (2.10)  reduces  to  the  total  probability  of  errors  (wrong  decisions): 


R(f)  =  Pr{H1,H2}  +  Pr{H2,H1} 


Accordingly,  the  LRT  decision  rule  in  Equation  (2.11)  becomes 


(2.12) 


or  equivalently, 


p(X\H2)  a  p, 

P(X\H,)  |  P2 


(2.13) 


P(H2 \X)  |  P(H2\X)  (2.14) 

Hi 

which  is  obtained  by  substituting  into  Equation  (2.13)  a  posteriori  probability  equali¬ 
ties  P(Hi\X)  =  Fl  and  P(H2\X)  =  ~ Pi-  ■  Therefore,  to  minimize  the  total 

error  probability,  the  hypothesis  decision  is  made  simply  by  choosing  the  maximum 
a  posteriori  probability.  The  LRT  under  this  circumstance  becomes  the  “Maximum 
A  Posteriori  (MAP)”  test  [21]. 

Under  the  specialized  Bayes  decision  rule  in  Equation  (2.13),  the  observation  space 
is  bisected  into  two  complementary  subspaces:  Li  and  L2.  When  X  falls  into  Li,  Hi 
is  declared;  otherwise  H2  is  declared.  Decisions  can  be  wrong  in  two  possible  ways: 
X  G  Hi  but  it  falls  into  L2,  or  X  €  H2  but  it  falls  into  La.  The  Bayes  error  is  defined 
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as  the  sum  of  probabilities  of  the  above  two  decision  errors. 


=  Pi  [  p(X\Ih)  dX  +  P2  [  p(X\H2 )  dX 

J  L2  J  L/ 1 


(2-15) 


—  P\€i  +  P2e  2 


It  can  be  proved  [3]  that  the  Bayes  error  e  is  the  lowest  achievable  total  probability 
of  error.  It  can  also  be  shown  [3]  that  the  Bayes  error  e  is  bounded  by  the  so-called 
“Chernoff  bound”  eu: 

-  p;pf  I {p(X\Hi))s  (p(X\H2))'- iX  VO  <  s  <  1 

(2.16) 

where  s  =  soptimUm  minimizes  eu. 

When  bothp(X|i/1)  and p(X\H2)  are  Gaussian,  the  integration  in  Equation  (2.16) 
will  lead  to  a  closed-form  expression.  Suppose  under  Hi  and  H2,  X  obeys  N(Mi,  Ei) 
and  N(M2,  E2),  respectively  (hereafter  N(M,  E)  denotes  a  Gaussian  distribution  with 
the  mean  vector  M  and  the  covariance  matrix  E).  Then  Equation  (2.16)  reduces  to 


eu{s)  =  P*  P]~s  e~M 


(2.17) 


where 


Ms)  =  “  M'^  +  (!  -  s)S2 r\M2  -  Mi)  + 

z  2  |Ei|s  IE2I1  s 

(2.18) 

is  called  the  “Chernoff  distance” . 

If  we  loosen  the  selection  criterion  of  s,  and  just  let  s  =  1/2,  the  Chernoff  bound 
specializes  to 
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£,(1/2)  =  V^-n/2) 


(2.19) 


where 


1  y  _i_  y1  i  I  Si+s2 1 

Ml/2)  =  -AM2  -  m1)t[-4-^]-'(M2  -  M,)  +  -In  2  ' 

o  z  2  y  |Si|  |E2| 

(2.20) 

is  called  the  “Bhattacharyya  distance” .  Actually,  if  both  classes  obey  Gaussian  dis¬ 
tributions  and  their  covariance  matrices  are  identical,  i.e.,  =  E2  =  E,  it  can  be 

easily  shown  that  soptjm„m  is  indeed  1/2.  In  this  case,  the  Bhattacharyya  distance 
in  Equation  (2.20)  is  further  simplified  to 


/A  1/2)  =  \{M2  -  M1)tE~1(M2  -  Mr)  (2.21) 

which  corresponds  to  the  lowest  Chernoff  bound.  In  Subsection  2.3.3,  we  will  relate 
the  Bhattacharyya  distance  with  the  class  separability  criterion  we  use  in  the  thesis. 

2.3.3  Feature  Extraction  and  Linear  Classifier 

Feature  extraction  is  to  choose  those  components  that  are  most  effective  for  separating 
classes.  It  is  a  process  of  converting  the  original  vector  Y  into  a  lower-dimensional 
feature  vector  Z.  Then  Z,  rather  than  Y,  is  fed  into  the  hypothesis  decision  rule. 
The  selection  of  Z  is  crucial  to  the  classifier  design:  if  Z  shows  significant  difference 
from  one  class  to  another,  the  classifier  has  a  good  performance.  The  selection  rule 
is  based  on  a  class  separability  criterion. 

Theoretically  speaking,  the  Bayes  error  e  as  in  Equation  (2.15)  is  the  ideal  crite¬ 
rion  for  class  separability,  but  a  major  disadvantage  is  that  an  explicit  mathematical 
expression  of  the  Bayes  error  is  not  available  except  for  a  very  few  special  cases  [3]. 
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Even  for  Gaussian  distributions,  calculation  of  the  Bayes  error  requires  a  numerical 
integration  unless  the  covariance  matrices  are  identical.  Thus  we  cannot  expect  in¬ 
sightful  theoretical  development  by  using  the  Bayes  error  as  the  feature  extraction 
criterion.  Corresponding  to  the  Bayes  error  as  the  ideal  criterion,  a  posteriori  prob¬ 
abilities  as  in  Equation  (2.14)  are  the  ideal  features.  Unfortunately,  in  practice,  a 
posteriori  probabilities  are  hard  to  obtain,  and  their  estimates  often  have  severe  bi¬ 
ases  and  variances  [3].  Thus  feature  extraction  based  on  the  Bayes  error  criterion 
tends  not  to  be  a  practical  solution. 

We  need  simpler  criteria  based  on  clear  physical  notions  and  associated  with  sys¬ 
tematic  feature  extraction  algorithms.  The  process  of  establishing  class  separability 
criteria  and  using  them  to  extract  features  is  called  “discriminant  analysis”  [3].  In 
the  thesis,  we  use  a  criterion  that  is  based  on  a  function  of  scatter  matrices  [3].  This 
criterion  is  formulated  by  a  within-class  scatter  matrix  and  a  between-class  scatter 
matrix.  The  within-class  scatter  matrix  AW_Y  depicts  the  scatter  of  samples  around 
their  respective  class  means: 

2  2 

A,.r  =  E  aE[(V  -  Mi)(Y  -  M,)t\H,]  =  (2.22) 

2=1  i=  1 

where  Y  is  the  data  column  vector;  Pi  =  is  the  a  priori  probability  of  class 

*5  Mi  =  E[Y\Hi ]  is  the  mean  vector  of  Y  in  class  i]  E*  is  the  covariance  matrix  of  Y 
in  class  i. 

The  between-class  scatter  matrix  Ab_Y  measures  the  “distance”  between  the  two 
classes: 


Ab_y  =  -  Mo ){Mi  -  Mof] 

i=i  (2.23) 

=  PiP2(M2  —  M1){M2  —  M\)T 


where  the  overall  mean  M0  =  ZlAPiMt) 
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The  class  separability  metric  is  defined  as 


Jy  =  tr{A-\YAb_Y )  =  PXP2  tr([J2(Pi^i)]-1(M2  -  M1)(M2  -  M if) 

1=1  (2.24) 

Trace  of  a  matrix  equals  the  sum  of  the  eigenvalues.  So  Equation  (2.24)  implies 
that  good  separability  requires  a  large  between-class  scatter  and  a  small  within-class 
scatter. 

It  is  worthwhile  to  relate  Jy  to  the  Bhattacharyya  distance  in  Equation  (2.20) 
under  a  special  condition.  When  Y  obeys  Gaussian  distribution  with  £i  =  £2  =  £ 
and  Pi  —  P2  —  1/2,  Jy  is  simplified  to 


Jy  =  \tr(E~1{M2  -  Mi)(M2  -  Mif)  =  \{M2  -  Mx)TYr\M2  -  Mx) 

ft 

(2.25) 

Comparing  Equation  (2.25)  with  Equation  (2.21),  we  see  that  under  the  conditions 
of  Gaussian  distributions,  equal  covariance  matrices,  and  equal  a  priori  probabilities, 
the  class  separability  Jy  and  the  Bhattacharyya  distance  fi  ( 1/2)  are  the  same  except 
for  a  constant  coefficient. 

With  the  class  separability  set  up,  let  us  now  establish  the  feature  extraction 
algorithm.  We  consider  a  linear  mapping  Z  =  CTY  that  transforms  the  observation 
vector  Y  into  a  lower-dimensional  feature  vector  Z: 


Z  -  CtY  (2.26) 

where  VisnxljZismxl;  the  transformation  matrix  C  is  n  x  m,  with  m  <  n 
for  dimension  reduction.  We  are  justifiably  concerned  about  any  possible  loss  of  class 
separability  information  induced  by  dimension  reduction,  but  it  will  be  shown  as 
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follows  that  under  the  separability  criterion  defined  by  Equation  (2.24),  the  dimension 
of  Z  is  lowered  to  one  without  inducing  any  loss. 

The  class  separability  metric  in  the  Z-space  is 

Jz  =  tr(A~1_zAb_z)  (2.27) 

where  by  the  definitions  of  within-class  and  between-class  matrices  as  in  Equation  (2.22) 
and  Equation  (2.23),  we  have 

A  w_z  =  CtAw_yC  (2.28) 

Ab_z  =  CTAb_YC  (2.29) 

Now  let  us  find  the  optimum  transformation  matrix  Coptimum  such  that  Jz  achieves 
its  maximum.  Without  considering  the  specific  structure  of  matrix  Ab_Y,  the  solution 
of  C  is 


^optimum  —  [bVl:by2*  *  '  '  ;bym]  (2.30) 

where  VY1,  VY2, ...  ,  VYm  are  the  eigenvectors  of  matrix  A^yA^  associated  with  its 
largest  m  non-zero  eigenvalues:  Ai,  A2,  ,  \m  in  the  descending  order.  JY  can 

be  expressed  as  the  sum  of  all  eigenvalues,  while  Jz  as  the  sum  of  those  largest  m 
eigenvalues: 


71 


i— 1 

(2.31) 

m 

y ^ 

!=1 

(2.32) 
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Taken  at  one  extreme,  when  m  =  n  and  C  is  an  invertible  matrix,  it  can  be  proved  [3] 
that 


Jz  =  JY  (2.33) 

i.e.,  a  non-singular  transformation  always  preserves  class  separability. 

By  its  definition  in  Equation  (2.23),  however,  Ab_y  is  special  in  that  its  rank 
equals  one.  Accordingly,  the  rank  of  A~1_YAb_y  also  equals  one,  so  the  matrix  has 
only  one  non-zero  eigenvalue.  Therefore  m  —  1,  and  Copt,mum  reduces  to  a  column 
vector,  which  is  the  eigenvector  Vy i  [3]  of  A^1_YAb_y  associated  with  its  only  non-zero 
eigenvalue  Ai: 


C optimum  —  Vy\  —  (3Aw]y(M2  —  Mx)  (2-34) 

where  /?  is  an  arbitrary  non-zero  constant  coefficient. 

It  is  thus  shown  that  under  the  separability  criterion  in  Equation  (2.24),  all  of  the 
class  separability  information  is  contained  in  one  eigenvector.  The  original  observation 
vector  Y  is  transformed  into  a  scalar  feature  z.  Nonetheless,  we  still  have  Jz  =  Jy, 
i.e.,  no  separability  information  is  lost.  The  scalar  feature  z  needed  for  classifying  Y 
is 


*  =  =  V&V  =  P(M2  -  M^A-iyY  (2.35) 

noting  that  Aw_y  is  symmetric.  Equation  (2.35)  represents  the  sufficient  statistic  of 
a  linear  classifier.  Once  the  statistic  is  chosen,  the  threshold  can  be  determined  by 
minimizing  the  total  cost  (in  the  Bayesian  sense),  or  by  satisfying  a  prescribed  false 
alarm  probability  (in  the  Neyman-Pearson  sense)  [21],  [22],  [23]. 

It  is  noted  that  the  Bayes  likelihood  ratio  test  is  always  the  optimum  classifier 
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since  it  minimizes  the  cost  or  probability  of  error.  To  construct  the  likelihood  ratio, 
one  has  to  estimate  the  conditional  probability  density  function  for  each  class  using  a 
finite  number  of  samples.  The  estimation  procedure  is  generally  very  complicated  and 
requires  a  large  number  of  samples  to  give  accurate  results  [3].  Even  if  the  conditional 
densities  can  be  estimated,  the  likelihood  ratio  test  is  often  difficult  to  implement: 
time  and  storage  requirements  may  be  excessive  [3].  We  are  therefore  often  led  to 
consider  a  simpler  procedure  for  designing  a  classifier.  In  particular,  we  may  specify 
the  mathematical  form  of  the  classifier,  leaving  only  a  finite  set  of  parameters  to  be 
determined. 

The  most  common  choices  are  linear  or  quadratic  classifiers.  The  linear  classifier 
discussed  above  becomes  the  Bayes  classifier  (i.e.,  the  optimum  classifier)  only  for 
Gaussian  distributions  with  equal  covariance  matrices  [3].  When  those  conditions  are 
not  met,  the  linear  classifier’s  performance  will  be  inferior  to  that  of  the  Bayes  classi¬ 
fier.  In  practice,  however,  the  linear  classifier’s  simplicity  and  robustness  compensate 
for  its  loss  in  performance  [3].  In  the  thesis,  we  adopt  the  class  separability  metrics 
Jy  and  Jz  (also  called  the  “Fisher  criterion”  [3])  as  just  reviewed.  They  represent 
clear  physical  notions  and  lead  to  a  linear  feature  extraction  algorithm. 
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Chapter  3 


Mingled  Spectrum  Principle 


3.1  Mingled  Spectrum  Recorded  by  a 
Moving  Platform 


3.1.1  Formula  Derivation 


Ocean  process  X(t,r) 


(to,  To) 


Y(t) 


speed=U 


(to+T,  To+TU) 


Figure  3.1:  A  line  AUV  survey. 


Oceanographic  processes  vary  both  in  time  and  space.  Under  the  assumptions  of  tem¬ 
poral  stationarity  and  spatial  homogeneity,  an  oceanographic  field  can  be  described 
by  its  frequency-wavenumber  spectrum.  When  an  AUV  (or  some  other  moving  plat¬ 
form)  makes  a  survey  in  the  field,  it  records  a  time  series  of  the  measured  quantities, 
e.g.  flow  velocity.  The  time  series  mixes  temporal  and  spatial  variations  of  the  sur¬ 
veyed  field.  The  corresponding  spectrum  therefore  mingles  the  spectral  information 
in  time  and  space,  which  we  call  a  ’’mingled  spectrum”. 

In  this  chapter,  we  present  the  mingled  spectrum  formula.  Suppose  an  AUV 
conducts  a  line  cruise  at  a  speed  of  u,  as  shown  in  Figure  3.1.  Denote  the  field  on 
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the  survey  line  as  X(t,r),  and  the  time  series  recorded  by  the  AUV  as  Y(t).  The 
autocorrelation  function  of  Y ( t )  is 


RY(r)  =  E[Y(t0)-Y(t0  +  r)] 

=  E[X(t0,r0)  ■  X(t0  +  t,  r0  +  ur)]  (3.1) 

=  Rx(r,  ut ) 

Then  the  Power  Spectrum  Density  (PSD)  of  Y(t),  i.e.,  the  mingled  spectrum,  is 
the  Fourier  transform  of  Ry(t )  (by  the  Wiener-Khinchine  theorem  [24]): 


/CO 

RY{r)e-j2^TdT 

-oo 

/CO 

Rx  (r,  ur)e~j2nfT  dr 

•OO 


(3.2) 


For  the  temporal-spatial  process  X,  its  autocorrelation  function  Rx(t,ut)  and 
its  PSD  Sx(r],  u)  are  a  Fourier  transform  pair  (also  by  the  Wiener-Khinchine  theo¬ 
rem  [24]): 


/oo  roo 

/  Sx{r}^)ej2^Te-j27U/UTdr]du  (3.3) 

oo  J  —  OO 

where  77  is  the  temporal  frequency,  and  v  =  k/(2Tr)  is  the  spatial  frequency.  Note 
that  A;  is  a  one-dimensional  wavenumber  in  the  direction  of  AUV’s  line  survey. 
Incorporating  Equation  (3.3)  into  Equation  (3.2),  we  have 
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/oo  poo  poo 

/  /  SxCv,  u)e~^e-^UTeJ^dTdridv 

-oo  J  —oo  J  —  oo 

/OO  /'OO  /'OO 

/  Sx(r},v)dr}dv  /  g-j'2iTT(/+Mt-jj)^r 

-oo  J —oo  J —oo 

/OO  /'OO 

/  Sx(?7,  ^)<5(/  +  ^  -  v)dri 

-OO  J  —OO 

/OO 

5x((/  +  H>^)^ 

-OO 


Figure  3. 2:, Illustration  of  the  derivation  of  Sy  from  Sx- 


Hence  the  mingled  spectrum  SV(/)  is  the  integration  over  v  of  Sx(v ,  on  a  line 
defined  by  rj  =  f  +  uu,  as  illustrated  in  Figure  3.2.  The  integration  line’s  slope  equals 
the  reciprocal  of  AUV  speed  u.  The  integration  line’s  intercept  on  the  77-axis  equals 
/.  Equation  SY{f) \u  -  Sx((/  +  w),  v)dv  thus  concisely  reveals  the  relationship 

between  the  “AUV-seen”  mingled  spectrum  SY{f)  and  the  original  temporal-spatial 
spectrum  Sx(v>  v). 

It  is  worthwhile  to  clarify  the  connection  between  the  “line”  PSD  Sx{ ri,  v)  used 
above  and  the  “complete”  PSD  SXj3d(v,  ^2)  or  its  variants.  The  discussion  is 
in  Appendix  A. 
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3.1.2  Special  Cases 

1.  Platform  speed  u  =  0. 


When  u  =  0,  Equation  (3.4)  is  simplified  to 


SY(f)\u= o=  Sx(f,  v)dv  (3.5) 

J  —  oo 

Referring  back  to  Figure  3.2,  u  =  0  makes  the  integration  line  perpendicular 
to  the  77-axis.  Since  the  measurement  platform  is  spatially  fixed,  the  mingled 
spectrum  at  any  frequency  /  is  obtained  by  integrating  the  temporal-spatial 
spectrum  Sx  over  the  spatial  frequency  u  at  77  =  /.  This  case  applies  to  a 
mooring  when  not  considering  any  advective  ambient  current. 

2.  Field  is  temporally  frozen. 

When  the  field  experiences  no  temporal  variation,  its  temporal-spatial  spectrum 
Sx  becomes  an  “impulse  fence”  [25]  on  the  i/-axis  (as  shown  in  Figure  3.3): 


Sx{rj,v)  =  Sxo(u)5(r]) 


(3.6) 


where  6(77)  is  the  unit  impulse  function  [26].  Then  Equation  (3.4)  becomes 


/OO 

Sxo(v)5(f  +  vu)di 

OO 


(3.7) 


By  variable  replacement  f;  ~  vu,  we  have 


Sy(/)l«  =  -  r  SM-m  +  f)di 
«  J-00  u 


iSxo(--) 

u  u 


(3.8) 


42 


Figure  3.3:  Derivation  of  Sy  from  Sx  for  a  temporally  frozen  field. 


The  mingled  spectrum  derivation  in  this  special  case  is  illustrated  in  Figure  3.3. 
Note  that  when  process  X  is  real,  Sxo(v)  is  an  even  function:  Sxo(—£)  — 
Sx o(f)-  Then  Equation  (3.8)  can  be  written  as 


Sk(/)I„  =  N*  o(-)  (3.9) 

u  u 

Equation  (3.9)  is  actually  an  expression  corresponding  to  the  Taylor’s  hypoth¬ 
esis  [27],  [9],  [28],  [29].  The  Taylor’s  hypothesis  is  widely  used  in  experimental 
turbulence  studies  [30],  [31],  [32],  [33],  [34],  and  also  in  ocean  acoustics.  The  hy¬ 
pothesis  states  that  if  the  traversing  speed  u  of  the  probe  is  large  enough,  we  can 
assume  that  the  turbulence  field  is  “frozen”,  i.e.,  no  temporal  variation  occurs 
within  the  measurement  duration.  Thus  although  the  probe  actually  records 
a  time  series  p(t),  it  can  be  transformed  to  a  spatial  series  q(x )  by  replacing 
t  by  |  [9].  The  Taylor’s  hypothesis  is  also  known  as  the  frozen-turbulence 
approximation  [28]. 

3.  Applied  to  a  nondispersive  plane  wave. 

As  an  anisotropic  process,  a  nondispersive  plane  wave’s  temporal-spatial  spec¬ 
trum  Sx  is  a  skewed  “impulse  fence”  [25]  (as  shown  in  Figure  3.4): 
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Sx(v,  v)  =  SXi{r])S{uc  -  77) 


(3.10) 


where  c  is  the  phase  speed  (also  group  speed  due  to  nondispersion)  of  the  plane 
wave,  and  it  is  assumed  herein  that  c  >  the  platform  speed  u  >  0.  Then  Equa¬ 
tion  (3.4)  becomes 


/OO 

Sx\ {f  +  vu)8(vc  -  (/  +  vu))dv 

■OO 


OO 

1  sxl(^L) 


(3.11) 


c  —  u 


c  —  u 


Figure  3.4:  Derivation  of  Sy  from  Sx  for  a  nondispersive  plane  wave  field. 

The  mingled  spectrum  derivation  in  this  special  case  is  illustrated  in  Figure  3.4. 
When  the  plane  wave  has  only  one  single  frequency  component  /0,  i.e.,  it  is  a 
tone  signal,  we  have 


Sxi(v)  =  a6(rj-  f0) 


(3.12) 


where  a  is  the  magnitude  coefficient.  Then  Equation  (3.11)  is  simplified  to 
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(3.13) 


Sy(f)\u  = 


-2-S(-SL 

c  —  u  c  —  u 


-fo) 


Letting  ^  —  /0  =  0,  we  find  the  frequency  of  the  signal  “seen”  by  the  moving 
platform: 


7/ 

/  =  (!--)/»  (3.14) 

This  is  actually  the  Doppler-shifted  frequency  recorded  by  a  moving  receiver  [35]. 

3.1.3  Differences  from  Doppler- Shifted  Spectrum  Method 

Comparing  the  presented  mingled  spectrum  principle  with  the  existing  work  of  Doppler- 
shifted  spectrum  as  introduced  in  Section  2.2,  we  observe  the  following  differences: 

1.  Different  requirements  on  isotropy. 

Doppler-shifted  spectrum :  Equation  (2.3)  assumes  horizontal  isotropy,  so  only 
the  radial  wavenumber  k  —  y/k\  +  k%,  rather  than  ki  together  with  k2,  is  needed 
for  the  spatial  axis.  For  anisotropic  fields,  therefore,  this  equation  cannot  be 
directly  applied. 

Mingled- spectrum:  Equation  (3.4)  does  not  require  the  field  to  be  isotropic.  It 
operates  on  one-dimensional  wavenumber  u  [v  =  |^).  This  equation  is  readily 
applicable  to  isotropic  or  anisotropic  fields. 

2.  Different  formats  and  conciseness. 

Doppler-shifted  spectrum:  the  integration  is  over  three  regions,  bounded  by  four 
lines.  Factors  [ k 2  —  (5£^;)2]~^  and  [ k 2  —  (5i^£)2]_^  in  the  integrands  further 
complicate  inspection  and  computation. 
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Mingled- spectrum,  the  integration  is  constrained  by  one  single  line.  The  inte¬ 
grand  is  simply  the  frequency-wavenumber  spectrum,  with  no  additional  factor. 
The  conciseness  facilitates  computation,  inspection,  and  interpretation. 

3.  Different  perspectives. 

Doppler-shifted  spectrum :  developed  from  the  perspective  of  the  Doppler  effect 
on  elementary  waves. 

Mingled- spectrum,  developed  from  the  perspective  of  temporal-spatial  autocor¬ 
relation  functions  and  by  utilizing  the  Wiener-Khinchine  theorem. 

As  detailed  in  Subsection  3.1.1,  our  mingled  spectrum  formula  (  Equation  (3.4))  is 
developed  directly  from  a  line  survey  in  a  temporal-spatial  field.  Consequently,  the  re¬ 
sultant  formula  is  not  constrained  by  isotropy  requirement.  Furthermore,  the  formula 
is  concise  and  embodies  a  clear  physical  interpretation.  Versatility  for  anisotropy,  ease 
for  inspection,  and  simplicity  of  computation  add  to  the  usefulness  of  the  presented 
formula. 


3.2  Utilization  for  AUV-Based  Classification 

Let  us  first  look  at  two  simple  fictitious  temporal-spatial  fields  for  the  purpose  of 
demonstration.  Their  rj-u  spectra  are  expressed  in  Equation  (3.15)  and  Equation  (3.16), 
and  displayed  in  Figure  3.5.  The  rj-u  spectrum  of  field  No.  2  is  just  a  transposition 
of  that  of  field  No.  1.  In  both  spectra,  the  range  of  frequency  rj  is  -1  Hz  ~  1  Hz  while 
the  range  of  wavenumber  v  is  -1  m-1  ~  1  m-1. 


Sx\{n,v) 


‘2,ko7j\gu\ 


-e  T?1 


(3.15) 


Sx2  (v,  v) 


2  7T (7^20^2 


(3.16) 
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where  r)0  =  u0  =  0.5 j  av i  =  =  0.2,  av 2  =  cx„i  =  0.1 


PSD  of  fictitious  field  No.  1 


-1  -0.8-0.6-0.4-0.2 


Temporal  frequency  (s  ) 


PSD  of  fictitious  field  No.  2 


-1  -0.8-0.6-0.4-0.2 


Temporal  frequency  (s  ) 


Figure  3.5:  Two  fictitious  rj-u  spectra. 


Suppose  an  AUV  is  to  classify  the  above  two  fields  based  on  a  line  survey  as 
illustrated  in  Figure  3.1.  It  is  noted  that  a  line  cruise  is  the  most  common  survey 
mode  for  an  AUV.  Its  advantage  is  long-distance  coverage,  thus  increases  chances  of 
finding  ocean  processes  of  interest.  During  such  a  survey,  the  vehicle  records  a  time 
series  of  the  concerned  quantity,  like  flow  velocity.  Based  on  distinct  spectra  of  the 
time  series  associated  with  different  processes,  the  AUV  carries  out  classification. 

, .  The  spectrum  of  the  AUV-recorded  time  series  is  just  the  mingled  spectrum  for¬ 
mulated  in  Subsection  3.1.1.  The  mingled  spectrum,  rather  than  the  field’s  original 
frequency-wavenumber  spectrum,  is  the  information  resource  for  spectral  classifica¬ 
tion,  because  time  and  space  are  mixed  in  the  AUV’s  record.  As  revealed  by  Equa¬ 
tion  (3.4)  and  Figure  3.2,  time-space  mixing  is  tuned  by  the  AUV  speed  u:  the 
integration  is  constrained  by  a  line  whose  slope  equals  l/u.  '■  ■  ■ 

For  the  two  fictitious  fields  whose  original  temporal-spatial  PSDs  are  given  by  Eq.ua- 
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Field  No.  1  Field  No.  2 


Figure  3.6:  Derivation  of  mingled  spectra  from  the  two  fictitious  rj-v  spectra. 

tion  (3.15)  and  Equation  (3.16),  their  mingled  spectra  under  some  AUV  speed  can  be 
easily  computed  by  Equation  (3.4),  as  illustrated  by  Figure  3.6.  At  a  series  of  vehicle 
speeds,  the  two  mingled  spectra  are  shown  in  Figure  3.7.  The  observation  is:  the 
two  mingled  spectra  may  appear  more  alike  or  more  distinct  depending  on  the  AUV’s 
cruise  speed.  Due  to  the  “transposition”  relation  between  the  two  hypothesized  77  —  u 
spectra,  their  mingled  spectra  are  identical  when  the  AUV  cruises  at  a  speed  of  1  m/s 
(the  third  panel  of  Figure  3.7).  This  would  obviously  prohibit  classification.  At  other 
speeds  of  0.5  m/s  (the  second  panel)  and  2  m/s  (the  fourth  panel),  however,  the  two 
processes  are  classifiable  since  their  mingled  spectra  show  difference.  A  quantitative 
metric  for  separability  will  be  given  in  Section  4.3. 

Our  goal  is  not  trying  to  reconstruct  the  field  [36],  [37],  [38]  or  its  original  spec- 
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temporal  frequency  (s-1) 

Figure  3.7:  Mingled  spectra  of  the  two  fictitious  fields. 


trum  [39],  but  to  classify  the  fields  by  the  difference  between  their  mingled  spectra 
acquired  by  an  AUV.  From  this  perspective,  we  are  to  utilize  the  mingling  of  time  and 
space  to  the  advantage  of  classification,  rather  than  regarding  the  Doppler  effect  as  a 
contaminating  factor  [40].  In  the  following  chapters,  we  will  develop  an  AUV-based 
spectral  classifier,  and  demonstrate  that  the  AUV  speed  can  be  optimized  to  high¬ 
light  distinct  features  of  processes.  We  will  use  two  real  examples  of  ocean  processes: 
ocean  convection  and  internal  waves. 


Chapter  4 


Spectral  Feature  Classification  of 
Different  Processes  Using  an  AUV 

4.1  Classifier  Architecture 

Now  we  apply  Equation  (2.35)  to  spectral  feature  classification  using  an  AUV.  In  the 
thesis,  we  study  two-class  problems.  The  classifier’s  architecture  is  illustrated  in  Fig¬ 
ure  4.1.  As  noted  in  Section  3.1,  The  AUV’s  measurement  Y(t)  mingles  temporal 
and  spatial  variations  of  the  field.  Hence  its  Power  Spectrum  Density  (PSD)  SY(f )  is 
a  mingled  spectrum.  It  is  related  to  the  field’s  temporal-spatial  PSD  Sx{i 7,  v)  by  the 
mingled  spectrum  formula  (Equation  (3.4)).  Classification  relies  on  the  “distance” 
(i.e.,  the  spectral  separability  to  be  given  in  Section  4.3)  between  the  two  mingled 
spectra  Syi(f)  and  SY2(f)-  The  AUV  speed  tunes  this  distance. 

From  AUV’s  measurement  Y(t),  we  obtain  the  estimate  of  its  PSD,  SY(f).  Here¬ 
after  we  denote  the  true  PSD  as  SY(f)  and  its  estimate  as  Sy(f )  (for  class  1  and  2, 
footnotes  “1”  and  “2”  are  added  for  distinction).  PSD  estimate  Sy(f )  is  the  input 
to  the  classifier.  We  use  the  Fourier  method  for  spectrum  estimation,  so  SY(f)  is 
given  at  a  series  of  discrete  frequencies.  Thus  the  PSD  estimate  is  expressed  as  a 
random  column  vector  SY(k),  k= 0,  1,  . . . ,  N  —  1,  where  N  is  the  total  number  of  fre¬ 
quency  points.  Then  for  classification,  the  scalar  feature  described  by  Equation  (2.35) 
becomes 
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Ocean  process  X(t,r) 


Figure  4.1:  Diagram  of  an  AUV-based  spectral  classifier. 


2  =  VtSy 


(4.1) 


where  V  is  the  transformation  vector  to  convert  a  random  vector  Sy  to  a  random 
variable  z.  The  same  as  in  Equation  (2.34),  V  is  expressed  as  (the  coefficient  /3  is 
dropped  since  it  does  not  affect  the  classifier’s  performance): 


V  =  A~lY(My2  -  My i) 


(4.2) 


where 


Myi  —  E[Sy\Hi]  i~  1,2 


(4.3) 


is  the  class  mean  and 


2 

Aw_y  —  YlPiXyi) 

i=l 


(4.4) 
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is  the  within-class  scatter  matrix.  is  the  prior  probability  of  class  i,  and  E Yi  is  the 
covariance  matrix  in  class  i: 


Eyj  —  E[(Sy  —  MYi)(SY  —  MYi)T\Hi\  2=1,2  ,  (4.5) 


Figure  4.2:  Mechanism  of  feature  projection  when  the  two-dimensional  E Yi  is  diagonal 
with  cryj(l)  =  <Jyj(2)  (modified  on  the  basis  of  Figure  4-7  and  Figure  10-1  of  [3]). 

To  help  explain  the  mechanism  of  feature  projection  as  expressed  by  Equation  (4.2), 
let  us  consider  a  very  simple  case.  Suppose  vector  SY  has  only  two  components.  If 
furthermore,  EYi  is  diagonal  with  equal  elements  (aYi(  1)  =  <7^(2)),  AW_Y  will  be 
a  scaled  identity  matrix  as  expressed  in  Equation  (4.4).  Then  according  to  Equa¬ 
tion  (4.2),  the  transformation  vector  V  will  coincide  with  the  vector  of  difference 
between  mean  vectors:  My 2  —  MYi,  except  for  a  constant  coefficient  not  affecting  the 
classifier’s  performance.  Feature  projection  is  then  illustrated  in  Figure  4.2:  when 
the  clusters  of  vector  SY  are  projected  onto  V  =  MY2  -  Myi,  the  distance  between 
the  clusters  of  scalar  feature  z  is  maximized. 

Now  let  us  take  one  step  further.  If  Eyj  is  still  diagonal  but  with  unequal  elements: 
e-S->  aYi(  1)  >  cryi(2),  Aw  Y  will  no  longer  be  a  scaled  identity  matrix,  thus  will  play 
a  role  in  the  formation  of  V,  as  shown  in  Equation  (4.2).  As  illustrated  in  Figure  4.3, 
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When  for  both  classes,  (7y(1)>(Jy(2): 


Figure  4.3:  Mechanism  of  feature  projection  when  the  two-dimensional  E Yj  is  diagonal 
with  aYi(l)  >  crYi( 2). 

the  role  of  Am_y  is  to  rotate  V  such  that  V’s  alignment  tilts  away  from  the  S'y(l)- 
axis  but  towards  the  SY( 2)-axis.  This  tilting  represents  a  penalty  on  the  <Sy(l)-axis 
projection  because  the  uncertainty  of  SY(  1)  is  larger  than  that  of  SY (2).  The  distance 
between  the  feature  z  clusters  is  maximized  using  the  tilted  V. 

Finally,  the  scalar  feature  z  is  compared  with  a  threshold  to  make  the  classification 
decision: 


*  =  VtSy  ^  7  (4.6) 

Hi 

where  7  is  the  threshold.  As  will  be  discussed  in  Section  4.4,  the  threshold  is  deter¬ 
mined  by  minimizing  the  total  cost  or  probability  of  error  (the  Bayesian  criterion), 
or  by  satisfying  some  prescribed  false  alarm  probability  (the  Neyman-Pearson  crite¬ 
rion)  [21]. 
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4.2  Statistics  of  PSD  Estimate 


As  shown  in  Equation  (4.1),  the  transformation  vector  V  is  the  key  to  the  classifier. 
V  is  determined  by  the  statistics  of  Sy,  as  shown  in  Equation  (4.2).  Hence  we  first 
need  to  know  the  statistics  of  the  PSD  estimate  SY. 

4.2.1  PSD  Estimation  Method 

We  use  the  periodogram  to  obtain  the  PSD  estimate  SY.  The  estimation  is  carried 
out  in  the  following  steps: 

1.  Since  the  recorded  data  yw(t )  is  of  finite  duration  ([-f ,  f  ]),  it  is  equivalent  to 
a  windowed  section  of  the  original  process  realization  y(t): 


yw(t )  =  y{t)w(t ) 


(4.7) 


where  w{t)  is  a  boxcar  window,  being  unity  over  [~f ,  f]  and  zero  elsewhere. 


2.  Apply  the  Fourier  transform  to  yw(t ): 


Yw(f\T)  =  T[yw(t )] 


The  periodogram  is  defined  as 


(4.8) 


SUf\T)  =  ±\Yw(f\T)\2  (4.9) 

3.  Do  frequency-domain  smoothing,  and  take  the  smoothed  periodogram  as  the 
PSD  estimate: 
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(4.10) 


Sy(f)=  E  KSw(f  +  ^\T) 

n=-Ns 

where  2NS  +  1  is  the  total  number  of  frequency  points  involved  in  the  smoothing, 
and  h(n)  is  the  weighting  function. 

It  can  be  proved  [41],  [42]  that  the  ensemble  average  of  Sw(f\T)  approaches  the 
true  PSD  Sy(/)  as  T  approaches  infinity1: 

limT^E[Sw{f\T)}  =  SYU)  (4-11) 

i.e.,  the  periodogram  Sw(f\T)  is  an  asymptotically  unbiased  estimate  of  Sy(f).  When 
the  window  w(t)J s  length  T  is  so  large  that  its  bandwidth  is  much  smaller  than  that 
of  Sy(f),  the  bias  of  the  periodogram  will  be  very  small.  However,  this  unsmoothed 
estimate  is  very  noisy: 


aSw(f\T) 

Sy(f ) 


1 


(4.12) 


i.e.,  the  periodogram’s  standard  deviation  is  as  large  as  its  asymptotic  mean.  The 
purpose  of  frequency-domain  smoothing  in  the  above  step  3  is  to  lower  the  estimate’s 
variance,  but  at  the  cost  of  frequency  smearing.  For  instance,  if  the  smoothing  weights 
in  Equation  (4.10)  are  uniform  with  hn  =  2Ar1+1 ,  the  standard  deviation  of  Sy(f)  will 
be  reduced  by  a  factor  of  \/2Ns  +  1: 


1In  practice,  a  periodogram  estimate  is  often  computed  using  a  single  realization  (e.g.,  a  time 
series)  of  the  random  process.  For  the  time-based  (versus  ensemble-based)  periodogram  to  be  mean¬ 
ingful  for  representing  the  statistical  variation,  the  random  process  should  be  ergodic  [42],  [43],  [24], 
[44]- 
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It  should  be  noted  that 


•  In  the  above  step  1,  a  selected  window  function  can  be  intentionally  applied 
on  top  of  the  boxcar.  Common  windows  include  Bartlett  (triangle),  Hanning 
(raised  cosine),  and  Hamming  (modified  raised  cosine).  Fourier  transforms  of 
the  latter  windows  have  much  lower  sidelobes  than  a  boxcar  window,  but  with 
wider  mainlobes  which  compromise  the  frequency  resolution.  A  lower  sidelobe  of 
a  window  function  means  less  interference  from  nearby  frequencies  for  spectral 
estimation  at  any  desired  frequency. 

•  Besides  frequency-domain  smoothing,  time-domain  segmenting  can  also  be  ap¬ 
plied  to  reduce  the  estimation  variance.  In  this  manner,  the  original  time  series 
is  first  divided  into  segments.  Then  a  periodogram  is  calculated  for  each  seg¬ 
ment.  With  each  periodogram,  frequency-domain  smoothing  can  be  conducted 
as  described  in  the  above  step  3.  Finally,  the  average  of  all  smoothed  peri- 
odograms  is  taken  as  the  PSD  estimate  SY(f). 

4.2.2  Statistics  of  PSD  Estimate  Sy(f) 

Based  on  its  formulation  in  Equation  (4.10),  the  expectation  and  the  covariance 
function  of  SY(f )  are  [41] 


Ns  noo 

£[&(/)]=  £  M  Sy(0\W{(-(f  +  ~))\2d(  (4.14) 

- - AT  *J  — 00  J- 


n—~Nl 


Ns  Ns 


~  roo 

Cov[SY(fi),SY(f2)}  =  ^  ^2  hmhn\  Sy(£)W*(£  —  (fi  +  —  ))W(£  —  (/2  +  —  ))d£|2 

m=—Ns  n=—Ns  J  1  1 

(4.15) 


s  n—  —  iys 


where  W (/)  is  the  Fourier  transform  of  the  window  function  w(t );  hn  is  the  weight 
for  frequency-domain  smoothing.  To  ensure  that  Sy(f)  is  an  asymptotically  unbiased 
estimation  of  Sy(f),  it  should  be  constrained  that  (f^  \W(f)\2df)  (Y^n=-N$  hn)  =  1- 
The  variance  of  Sy(f)  can  be  obtained  directly  from  Equation  (4.15): 


Ns  Ns 

Var[SY(f)}  =  E  E^ 

m=—Ns  n=~Ns 


I  J"  Sr(pw*(?  -  </ + ^))w{s  -  (/ + ~))d(  |2 

(4.16) 


When  Sy(f )  is  smooth  across  the  smoothing  bandwidth  2Ap-!1 ,  Sy(f )  can  be 
pulled  out  of  the  integration  in  Equation  (4.16),  such  that 


{Ns  Ns  /»oo 

E  E  w"*(5  -(/ 

m——Ns  n——Ns  J 


(sAfjl 

v*!t 


j))w(a-u  +  ^))<ui  I2} 

(4.17) 


where  is  used  to  represent  the  quantity  in  the  curly  braces,  and  uefj  is  called 
the  “effective  number  of  degrees  of  freedom” .  The  cause  for  this  terminology  is  that 
Sy(f )  is  shown  [41]  [42]  to  approximately  obey  a  x2  distribution  with  2 ueff  degrees 
of  freedom2,  i.e., 


2**,/^  ~  X2(2 Veff)  (4-18) 

Another  important  observation  from  Equation  (4.15)  is  that  the  PSD  estimates  are 
approximately  uncorrelated  between  frequencies  farther  apart  than  Bw  (Bw  denotes 
the  bandwidth  of  W(f)): 


2The  real  y2  distribution  with  2i/e//  degrees  of  freedom  is  derived  from  a  complex  %2  distribution 
with  i /eff  degrees  of  freedom.  The  variance  of  each  component  equals  |  in  the  former  distribution, 
and  equals  1  in  the  latter  [41]. 
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Cov[SYi{fi),SYi{f2)}  ~  0  when  |/2  -  fi\>  Bw  *  =  1,2 

(4.19) 

To  give  an  idea  of  how  i/e^  is  related  to  window  shape  and  frequency-domain 
smoothing,  let  us  look  at  an  example.  When  a  boxcar  window  is  applied  to  the  time 
series  and  an  X-point  smoothing  (with  uniform  weights)  is  done  in  the  frequency  do¬ 
main,  we  have  uefj  =  y/L.  According  to  Equation  (4.17),  the  variance  of  SY  decreases 
by  a  factor  of  So  the  more  smoothing,  the  lower  the  estimation  variance.  The 
cost  paid,  however,  is  smearing  of  the  spectrum  within  the  smoothing  bandwidth.  Be¬ 
cause  of  this  smoothing,  the  total  number  of  frequency  points  (within  some  prescribed 
frequency  range)  that  provide  uncorrelated  PSD  estimates  are  effectively  reduced. 

4.3  Feature  Extraction  and  Spectral  Separability 

4.3.1  Computation  of  Feature  Transformation  Vector  V 

As  shown  by  Equation  (4.1),  the  transformation  vector  V  is  the  key  to  classification. 
To  compute  V  as  in  Equation  (4.2),  we  need  class  mean  vectors  My\  and  MY2,  as  well 
as  the  within-class  scatter  matrix  AW_Y.  Relating  Equation  (4.14)  to  Equation  (4.3), 
we  have 


AX 


Yl 


Ns  poo 

(!)  =  E[Sr{})m=  Y,  M  Syl(Or(f-(/+£))f 

- AT  J  ~  OO  J- 


dt 


n=—Nt 


(4.20) 


AX 


Y  2 


Ns  poo 

(/)  =  E[Sy(f)\H2]  =  Y  h-  SK2«)|W'(f-(/+£))|! 

- AT  J  ~  OO 


dt 


n=—N. 


(4.21) 


Matrix  AW_Y,  however,  is  composed  of  the  covariance  matrices  Eyi  and  Ey2,  as 
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shown  in  Equation  (4.4).  We  pick  frequency  points  with  an  interval  of  Bw,  so  the 
PSD  estimates  are  uncorrelated  according  to  Equation  (4.19).  Covariance  matrix  E Yi 
is  consequently  diagonal: 


E Yi  =  diag{Var[SYi(0)],  f/arfSy^l)],  •  •  •  ,  Var[SYi(N  -  1)]}  i  =  1,2 

(4.22) 

where  (0),  (1),  •  •  •  ,  (N  —  1)  denote  the  N  frequency  points.  By  Equation  (4.17)  (but 
replacing  Syi(f)  with  MYi(f)  to  take  into  account  the  window  effect),  the  diagonal 
elements  are  3: 


Var[SYi(f)]  = 


(■ MYi(f ))2 

veff 


i  =  1,2 


(4.23) 


With  Eyi  and  Sy2,  the  within-class  scatter  matrix  AW_Y  can  be  constructed 
by  Equation  (4.4)  (assuming  Px=  P2  =  |): 


Aw_y  =  y^(P;Eyj) 

i=l 

=  - - diag{(MYi( 0))  4- (My2(0))  ,  •••, 

tVeff 

(MY1(N-1))2  +  (MY2{N-1))2} 


Incorporating  Equation  (4.24)  into  Equation  (4.2),  we  have 


3The  zero-frequency  is  special:  Ear[5y(0)]  =  2- — P  '  .  This  specialty  is  properly  treated  in 
computations,  but  omitted  in  expressions  herein  for  the  sake  of  conciseness. 
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V  —  Aw1_y(MY2  —  My\) 

MY2(0)  —  Myx(  0) 

+  (My2(0))2]  (4.25) 

MY2{N  -  1)  -  Myi{N  -  1)  1T 
rfjy  [{MY1(N  -  l))2  +  ( MY2(N  -  l))2]  . 

4.3.2  Scalar  Feature’s  Statistics  and  Spectral  Separability 

Incorporating  Equation  (4.25)  into  Equation  (4.1),  we  have 


VTSy 


N- 1 

=  £ 


MY2(k)  —  Myi(h) 


k=0  2u, 


1 

'«// 


[(Myi(/c))  +  ( MY2(k ))  ] 


SY(k) 


Applying  Equation  (4.3),  we  obtain  the  mean  of  2: 


(4.26) 


mzi  =  VTE[SY\Hl] 
=  VTMYi 


N-l 


-  2  Veff 


MY2{k)  -  MY1{k) 
(Mvi(k))2+(MY2(k)) 


(4.27) 


;MYi(k )  *  =  1,2 


Applying  Equation  (4.22),  Equation  (4.23),  and  Equation  (4.25),  we  obtain  the 
variance  of  2: 


a^  =  VJEYiV 


.  r  MY2(k)  —  MYi(k)  -.2/,, 

~  V”’  S  (MnW)2+(My2(i))2^  ^  Yi(k) ) 


1,2 


(4.28) 


Accordingly,  the  2-domain  within-class  scatter  matrix  actually  reduces  to  a  scalar 
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(assuming  Px  =  P2  =  \): 


2 

<4._,  = 

i= i  (4.29) 

=  \{°2zl  +  <4) 

Now  let  us  examine  the  class  separability  metrics  in  the  Y  and  2-domains,  using 
the  general  definition  given  by  Equation  (2.24).  By  the  definition  in  Equation  (2.23), 
the  between-class  scatter  matrix  in  the  F-domain  is  (assuming  Pj  =  P2  =  5): 


Ab_y  —  ~(My2  ~  MYi){My2  —  Myi)T 


(4.30) 


and  the  between-class  scatter  matrix  in  the  2-domain  is  a  scalar: 


Ab_z  =  ^(mz2  ~  rrizi)2  (4.31) 

Applying  Equation  (2.24),  Equation  (4.24),  and  Equation  (4.30),  we  obtain  the 
E-domain  class  separability: 


JY  —  tr(Aw]yAb_Y ) 

=  ^  tr(A~1_Y(MY 2  —  MYi)(MY2  —  Myi)t) 

=  —  (MY2  —  MYi)t  Aw]y(MY2  —  My\ )  ^  ^ 

=  Veff  [MY2(k)  -  MY1(k)}2 

2  h  (MYi(k))2  +  {MY2(k))2 

Incorporating  Equation  (4.27)  and  Equation  (4.28)  into  Equation  (4.29)  and  Equa¬ 
tion  (4.31),  we  obtain  (noting  that  Aw ^  and  Ab_z  are  both  scalars)  the  2-domain  class 
separability: 
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(4.33) 


Jz  =  tr{A~\Ab^) 

_  (mz2  -  mzl)2 
2(^i  +  °l-i) 

_  Ueff  [My2(fc)  -  Myijk)]2 
2  £  {MY1(k))2  +  (MY2(k))2 

Comparing  Equation  (4.32)  and  Equation  (4.33),  we  see  that 


Jz  =  Jy  (4.34) 

This  verifies  the  general  rule  expressed  by  Equation  (2.33).  Under  the  definition 
given  by  Equation  (2.24),  separability  is  thus  preserved  after  dimension  squeezing 
from  an  TV  x  1  PSD  estimate  vector  Sy  to  a  scalar  feature  z. 

Let  us  make  some  observations  of  the  spectral  separability  metric  in  Equation  (4.33). 
The  numerator  indicates  the  difference  between  the  two  spectra.  Its  growth  en¬ 
larges  the  spectral  “distance”.  The  periodogram’s  variance  is  proportional  to  the 
square  of  the  spectrum  height,  as  shown  in  Equation  (4.17).  Thus  the  denomina¬ 
tor  in  Equation  (4.33)  actually  represents  the  spectrum  estimate’s  uncertainty.  Its 
growth  effectively  decreases  the  separability  because  we  have  less  trust  in  the  spectral 
difference.  Separability  increases  with  ueff,  but  we  should  keep  in  mind  the  involved 
cost.  If  uejf  >  1  is  realized  by  temporal  segmentation,  it  would  require  a  longer  data 
record.  With  longer  data,  however,  the  classifier’s  underlying  assumptions  of  tempo¬ 
ral  stationarity  and  spatial  homogeneity  would  become  more  challenged.  If  vejf  >  1 
is  realized  by  frequency-domain  smoothing,  it  would  smear  the  spectrum  over  the 
smoothing  bandwidth.  Within  some  prescribed  frequency  range,  the  smearing  would 
render  fewer  frequency  points  that  can  provide  uncorrelated  spectrum  estimation. 
Fewer  frequency  points  implies  lower  N  for  summation  in  Equation  (4.33),  which 
would  adversely  affect  separability. 

Another  note  is  that  so  far,  we  have  been  considering  the  periodogram’s  inherent 
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uncertainty  as  the  only  source  of  Var[SY}.  This  leads  to  a  very  simple  format  of  the 
denominator  in  Equation  (4.33).  When  external  uncertainties  are  added  in,  V  ar[SY) 
will  grow,  and  its  form  will  become  more  complicated.  The  necessity  and  methods  to 
treat  additional  uncertainties  will  be  presented  in  Chapter  7. 

4.4  Decision  Threshold  Setting  and  Classifier  Per¬ 
formance 

The  classification  decision  is  made  by  comparing  the  scalar  feature  with  a  threshold: 
h2 

z  \  7.  Once  the  sufficient  statistic  2  is  given  as  in  Equation  (4.1),  the  classifier’s 
Hi 

Pd-Pf  relationship  will  be  fixed  (Pd  denotes  the  detection  probability  and  PF 
denotes  the  false  alarm  probability) .  Here  we  have  adopted  the  terminology  in  signal 
detection  for  our  classification  scenario.  PF  is  defined  as  the  probability  of  deciding 
class  2  when  class  1  is  true,  i.e.,  PF  =  Pr{H2\Hi}.  PD  is  defined  as  the  probability  of 
deciding  class  2  when  class  2  is  indeed  true,  i.e.,  Pd  =  Pr{H2\H2}.  The  threshold  7  is 
to  be  determined  based  on  the  chosen  performance  criterion.  One  common  criterion 
is  to  minimize  the  Bayes  cost  [21].  Another  criterion  is  to  satisfy  some  prescribed 
Pf  [21], [22], 

According  to  Equation  (4.18),  the  PSD  estimate  SY  obeys  a  x2  distribution.  When 
uef  j  is  large,  we  know  that  a  x2  distribution  approaches  a  Gaussian  distribution,  as 
displayed  in  Figure  4.4.  For  the  sake  of  demonstration,  we  will  take  this  approximation 
in  the  remainder  of  this  section  (when  uejj  is  small,  the  concept  to  be  conveyed  is 
still  valid  but  closed- form  probability  calculations  are  no  longer  available).  Since 
z  =  VTSY  is  a  linear  transform  of  SY,  z  will  also  approximately  obey  a  Gaussian 
distribution: 


z\Hi~  N(mzi,o2zi)  i  —  1,2  (4.35) 


where  mzi  and  o2zi  are  z’s  mean  and  variance  in  class  i,  respectively.  Then  z’s  proba- 
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Figure  4.4:  A  x2  distribution  approaches  a  Gaussian  distribution  with  an  increasing 
number  of  degrees  of  freedom. 

bility  density  function  (PDF)  is  expressed  by 


Pz(z\Hi)  =  -=L—  exp(-^~  mzi^  )  i  =  1, 2  (4.36) 

y/2naZi  2a2zi 

Suppose  for  the  two  Gaussian  distributions,  mz i  =  -0.3,  mz2  =  0.6,  azl  =  0.2, 
and  cr^2  =  0.1.  The  two  PDFs  and  threshold  7  are  illustrated  in  Figure  4.5.  PF  is 
(z\Hx)  s  tail  probability  to  the  right  of  7,  while  1  —  P b  is  (2|FT2)’s  tail  probability  to 
the  left  of  7.  For  Gaussian  distributions,  we  have 


PF  =  0.5[1  -  signal)  erf( 7^)]  (4.37) 

PD  =  0.5[1  -  sign{^)  erf( j*)\  (4.38) 

where  7*  -  i  =  1, 2,  and  the  error  function  is  defined  as  erf(x )  =  ~  J0X  e-*2dt. 

As  7  varies,  P p  and  PF  vary  accordingly.  The  PF  —  PF  curve  is  also  referred  to  as 
the  Receiver  Operating  Characteristic  (ROC)  [21],  The  ROC  evaluates  a  classifier’s 
performance.  It  should  be  noted  that  the  ROC  is  determined  solely  by  the  PDFs 
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Figure  4.5:  Illustration  of  Gaussian  distributions  under  H\  and  H2,  and  definitions 
of  PF  and  Pp. 

Pz(z\Hi)  and  pz(z\H2),  not  by  the  threshold  7.  The  role  of  7  is  to  pick  one  particular 
spot  on  the  ROC  curve.  Generally,  a  smaller  7  corresponds  to  both  a  higher  PF 
and  a  higher  PD.  Choice  of  7  depends  on  the  user’s  preference  of  criterion.  Two 
common  criteria  are:  the  lowest  total  probability  of  error  (Bayesian  criterion),  and  a 
prescribed  PF  (Neyman- Pearson  criterion).  The  total  probability  of  error  is  defined 
a s  Pe  =  P\PF  +  P2(l  —  Pd)  where  Pi  and  P2  are  prior  probabilities  of  the  two  classes. 
In  the  sequel  we  assume  Pi  =  P2  =  |,  and  then  Pe  =  |(1  —  PD  +  PF ).  The  ideas  are 
demonstrated  by  the  ROC  curve  shown  in  Figure  4.6  that  corresponds  to  Figure  4.5. 
If  we  want  to  achieve  the  lowest  total  probability  of  error  P£,  7  should  be  set  to  0.17. 
On  the  other  hand,  if  we  prescribe  the  false  alarm  probability  as  PF  =  0.1,  then  7 
should  be  set  to  0.27  to  meet  the  requirement. 

Now  it  is  time  to  relate  the  classifier’s  ROC  to  the  class  separability  metric  Jy. 
According  to  Equation  (4.34),  Jz  =  Jy.  So  we  can  substitute  Jz  for  Jy.  Still  using 
the  two  PDFs  shown  in  Figure  4.5,  two  ROC  curves  corresponding  to  two  different 
values  of  Jz  are  compared  in  Figure  4.7.  The  curve  labeled  Jy  =  Jz  =  J  =  1.35  is  the 
same  as  in  Figure  4.6.  The  curve  below  has  Jy  =  Jz  =  J  =  0.675.  We  understand 
that  a  higher  ROC  curve  implies  a  better  classifier  performance,  no  matter  by  the 
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Figure  4.6:  Illustration  of  the  ROC  curve  and  the  relationship  between  the  threshold 
7  and  (PF,PD). 

Bayesian  criterion  or  the  Neyman-Pearson  criterion.  Thus  the  classifier’s  performance 
improves  with  class  separability  JY  —  Jz. 

4.5  Impact  of  AUV  Speed  on  Classifier  Perfor¬ 
mance 

As  pointed  out  in  Subsection  3.1.1,  when  an  AUV  makes  a  line  survey  in  field  X 
whose  temporal-spatial  PSD  is  Sx{r},v),  the  mingled  PSD  SY(f)  recorded  by  the 
AUV  is  dependent  on  the  vehicle’s  cruise  speed.  The  relationship  is  governed  by 
the  mingled  spectrum  formula  (Equation  (3.4)).  We  have  observed  in  Section  3.2 
that  the  mingled  spectra  of  two  different  ocean  processes  (SY i  v.s.  SY2)  may  appear 
more  alike  or  more  distinct  depending  on  the  AUV’s  cruise  speed.  The  quantitative 
metric  of  their  spectral  distance  is  expressed  by  JY  in  Equation  (4.32).  Note  that 
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Figure  4.7:  The  classifier’s  performance  improves  with  Jy  =  Jz. 

in  Equation  (4.32),  Mn  is  just  a  variant  of  Sn  when  taking  into  account  the  effect 
of  finite  data  length. 

As  shown  by  Figure  4.7,  the  linear  classifier’s  performance  relies  on  Jz  =  Jy.  We 
therefore  project  that  given  temporal-spatial  PSDs  of  two  different  ocean  processes, 
classification  can  be  better  conducted  at  some  AUV  speeds  than  others.  We  will 
next  introduce  two  kinds  of  ocean  processes,  and  study  the  AUV-based  classifier’s 
performance  as  a  function  of  the  vehicle  speed. 
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Chapter  5 


Ocean  Convection  and  Internal 
Waves 

5.1  Ocean  Convection 

5.1.1  Background 

Convection  is  the  transfer  of  heat  by  mass  motion  of  fluid  [35].  It  happens  when  the 
density  distribution  becomes  unstable  [8].  Open  ocean  convection  takes  place  at  only 
a  few  locations  around  the  world,  namely,  the  Labrador  Sea  [45],  [46],  the  Greenland 
Sea  [47],  [48],  Mediterranean  [16],  and  around  the  Antarctica  [49].  At  these  locations, 
strong  winter  cooling  of  the  surface  water  causes  it  to  become  denser  than  the  water 
beneath.  The  cooled  surface  water  sinks  and  mixes  with  deeper  water  which  enters 
the  global  ocean  circulation.  This  process  releases  heat  from  the  overturned  water  to 
the  atmosphere  and  thus  maintains  a  moderate  winter  climate  on  the  land.  Hence 
ocean  convection  is  an  important  mechanism  for  global  heat  transfer  [35],  [50],  [51]. 

At  those  open  ocean  convection  sites,  water  is  weakly  stratified  during  the  winter, 
i.e.,  there  exists  an  upper  mixed  layer.  When  winter  storms  set  in,  they  induce  intense 
heat  flux  from  the  sea  to  the  air.  Sea  surface  cooling  can  initiate  ocean  convection. 
In  convection  regimes,  the  water  column  overturns  in  numerous  convective  cells  (also 
called  plumes)  [52].  The  mixed  layer  thus  deepens,  and  gradually  the  convective 
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plumes  form  a  deep  “mixed  patch”  ranging  tens  of  kilometers  in  horizontal  scope  [52]. 
In  the  ending  stage,  with  the  cessation  of  winter  surface  heat  loss,  the  mixed  patch 
laterally  exchanges  water  with  its  surroundings.  The  mixed  water  thus  disperses 
under  the  influence  of  gravity  and  rotation,  leading  to  disintegration  of  the  mixed 
patch  and  re-occupation  by  stratified  water  [52], 

For  process  recognition,  a  very  useful  feature  of  convection  is  the  spatial  periodicity 
of  convective  cells  (as  will  be  shown  in  Figure  5.2).  Numerical  [53]  and  experimen¬ 
tal  [54]  studies  have  revealed  the  spatial  scale  of  convective  cells  as  a  function  of  the 
surface  heat  flux  and  the  Coriolis  frequency  (due  to  the  Earth’s  rotation). 

Consider  a  water  column  that  is  vertically  mixed  at  its  initial  state.  Heat  flux  of 
Q  is  applied  to  its  surface.  When  the  water  column  is  sufficiently  deep1 ,  the  evolving 
convective  layer  will  come  under  the  geostrophic  control  (due  to  the  Earth’s  rotation) 
as  time  approaches  /_1  where  /  is  the  Coriolis  frequency  [53].  Corresponding  to  the 
heat  flux  Q,  the  buoyancy  flux  B0  can  be  calculated  by 


Bo 


agQ 

pCp 


(5.1) 


where  a  =  2x  10~4  (°C)~1  is  the  water’s  thermal  expansion  coefficient;  g  =  9.81  m/s2 
is  the  acceleration  of  gravity;  p  is  the  water  density;  C/=3900  J/(kg  °C)  is  the  specific 
heat  capacity  of  water. 

Then  the  horizontal  scale  l  of  each  convective  cell  will  follow  [53]: 


(B-2) 

For  instance,  at  the  1998  Labrador  Sea  Experiment  site  (please  refer  to  Chapter  9), 
latitude  (j)  «  57°N.  Then  the  Coriolis  frequency  /  =  2f2sm(0)  ^  1.2  x  10-4  rad/s 
(fl  =  27r/day  is  the  angular  velocity  of  Earth’s  rotation).  Suppose2  the  heat  flux 

1i.e.,  depth  >  l  ~  (yj) 2 .  See  Equation  (5.2). 

2Surface  heat  flux  and  sea  water  density  values  cited  herein  are  measurements  during  AUV 
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Q— 300  W/m2  and  the  sea  water  density  p=1030  kg/m3,  then  the  buoyancy  flux 
B0  fa  1.5  x  10~7  m2/s3.  Accordingly,  a  convective  cell’s  horizontal  scale  is  (in  the 
calculation  by  Equation  (5.2),  the  value  of  /  should  be  in  unit  rad/s  [53]): 


l  ~  ( 


Bq _  /1.5  x  10  7  m2  s  3 

T3'  ~~  ^  (1.2  x  10-4  s-1)3 


1 

2 


290  m 


(5.3) 


Convective  cell’s  length  scale,  or  equivalently,  the  spatial  periodicity  of  convective 
cells,  is  what  the  AUV  utilizes  for  classification,  as  will  be  presented  in  Chapter  6. 


5.1.2  MIT  Ocean  Convection  Model 

Prof.  John  Marshall  and  his  group  at  the  MIT  Department  of  Earth,  Atmospheric, 
and  Planetary  Sciences  have  established  a  numerical  model  of  open  ocean  convec¬ 
tion  [53],  [55],  [56].  In  the  thesis,  we  use  this  model  to  find  the  temporal-spatial 
spectrum  of  convective  vertical  velocity. 

Heat  flux 

A 


Figure  5.1:  Illustration  of  the  convection  model  box.  Only  the  top  surface  is  subjected 
to  a  heat  flux.  The  dimension  is  200  x  200  x  35  with  a  grid  size  of  10  m. 

The  model  is  configured  as  a  box  as  shown  in  Figure  5.1.  Water  is  cooled  at  the 
top  surface.  There  is  no  normal  heat  flux  at  the  bottom  or  the  four  side  walls.  The 

Mission  B9804107  in  the  Labrador  Sea.  Please  see  Chapter  9. 
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mechanism  of  the  model  [57]  is  summarized  as  follows. 


•  Governing  equations. 

The  state  of  the  ocean  at  any  time  instant  is  characterized  by  the  three- 
dimensional  distribution  of  current  velocity  U  =  [u  v  w],  potential  temperature 
T,  salinity  S,  pressure  p  and  potential  density  p.  For  brevity,  the  potential 
temperature  and  the  potential  density  will  be  simply  referred  to  as  temperature 
and  density  in  the  following.  The  Boussinesq  approximation  [9]  is  applied  to 
the  model.  This  approximation  has  two  implications: 

1.  Fluid  is  incompressible.  This  leads  to  Equation  (5.4)  in  the  following. 

2.  Density  variation  is  negligible  in  momentum  equations  for  the  horizontal 
plane,  but  is  important  in  the  vertical.  Thus  the  density  variation  term  p' 
only  appears  in  Equation  (5.7),  but  not  in  Equation  (5.5)  or  Equation  (5.6). 

The  governing  equations  of  the  model  are  listed  as  follows, 


1.  Continuity  equation: 


du  dv  dw 
dx  dy  dz 


(5.4) 


where  u,  v,  and  w  are  the  two  horizontal  and  the  vertical  velocity  compo¬ 
nents. 

2.  Momentum  equations: 


Du  1  dp '  ,d2u  d2u .  d2u 

~Dt  +  ~ iv  =  l'hUx>- +  Thj^  +  ""a? 


(5.5) 


Dv  1  dp' 
Dt  +  po  dy 


,  ,d2v  d2v .  d2v 

+  fU  ~  Vh{W  +  Vvd* 


(5.6) 
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(5.7) 


Dw  1  dp'  p'  d2w  d2w  d2w 

Dt  p0  dz  p0  dx2  dy 2  dz 2 

where  po  is  the  static-state  reference  density;  p'  is  the  density  deviation 
from  p0;  g ^  is  the  buoyancy  forcing;  p'  is  the  pressure  deviation  from 
the  hydrostatic  pressure;  /  is  the  Coriolis  frequency;  vh  and  uv  are  the 
horizontal  and  vertical  eddy  viscosity,  respectively; 

3.  Heat  equation: 


DT  _  ^  (d2T  k  d2T ,  t  d2T 
Dt  ~  ^dx2  +  ~df  '  +  Kv Ik2 

where  T  is  the  temperature;  /c^  and  kv  are  the  horizontal  and  vertical 
diffusivity  of  temperature,  respectively. 

4.  Equation  of  state: 


P  —  Po[l  —  cx(T  ~  To)  +  p(S  —  aSo)]  (5.9) 

where  a  is  the  thermal  expansion  coefficient;  (3  is  the  saline  contraction  co¬ 
efficient.  To  and  So  are  the  static-state  reference  temperature  and  salinity, 
respectively. 

•  Boundary  conditions. 

The  boundary  conditions  of  the  model  are  as  follows, 

1.  Kinematic  boundary  conditions: 

-  At  top  and  bottom. 

Free  slip: 


du  dv 
dz  dz 


(5.10) 
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No  normal  flow: 


w  —  0 


(5.11) 


-  At  four  vertical  walls. 
Free  slip: 


du  dw 
dy  dy 


0  (for  x-directional  walls)  and 


dv  dw 
dx  dx 


0  (for  y-directional  walls) 


(5.12) 


No  normal  flow: 


v  =  0  (for  x-directional  walls)  and  u 


0  (for  y-directional  walls) 


(5.13) 


2.  Thermal  boundary  condition: 

-  At  bottom  and  four  side  walls. 
No  normal  heat  flux: 


(5.14) 


—  At  surface. 


k.cM  =  5-  (5.15) 

oz  p0 

where  Q  is  the  prescribed  surface  heat  flux. 

5.1.3  Model  Parameter  Setting 

In  the  thesis,  the  model  parameters  are  set  by  using  the  meteorological  and  hydro- 
graphic  data  acquired  in  the  1998  Labrador  Sea  Experiment  that  will  be  presented  in 
detail  in  Chapter  9.  In  this  way,  the  model  output  of  vertical  flow  velocity  can  be  rea¬ 
sonably  used  as  the  theoretical  template  for  matching  the  corresponding  measurement 
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made  by  the  AUV  (please  see  Chapter  9).  The  classifier  tests  will  be  presented  in  Sec¬ 
tion  9.5.  The  model  parameters  and  physical  parameters  are  listed  in  the  following 
two  tables,  and  explained  accordingly. 


Table  5.1:  Model  parameters 


Parameter 

Symbol 

Value 

Initial  reference  temperature 

T0 

Labrador  Sea  profile  * 

Initial  reference  salinity 

So 

Labrador  Sea  profile  * 

Surface  heat  flux 

Q 

300  W/m2  t 

Initial  velocities 

u0,  Vo,  Wo 

0 

Grid  length 

Ax,  Ay,  A z 

10  m  ° 

Number  of  grids  in  x  and  y 

NX,NV 

200  ° 

Number  of  grids  in  z 

Nz 

35  ° 

Internal  iteration  time  step 

At  Oration 

10  s 

Data  output  time  step 

AtouipUi 

30  s 

The  initial  temperature  and  salinity  profiles  are  assigned  values 
*  measured  by  the  AUV  during  Mission  B9804107  in  the  Labrador  Sea 
on  February  10,  1998. 


The  surface  heat  flux  is  assigned  a  value  obtained  from  the  measure- 
I  ment  on  board  R/V  Knorr  during  AUV  Mission  B9804107.  Details 
are  contained  in  Table  9.1  in  Section  9.2. 

Considering  the  characteristic  convective  wavelength  of  about  200  m  in 
the  350-m  mixed  layer  during  AUV  Mission  B9804107,  we  set  the  grid 
length  and  the  number  of  grids  such  that  there  can  be  held  about  10 
convective  cells  in  the  x  or  y  direction,  and  the  spatial  sampling  interval 
is  about  1/20  of  the  characteristic  wavelength.  Note  that  the  total 
number  of  grid  points,  Nx  x  Ny  x  Nz  is  constrained  by  the  processor 
capability.  The  job  defined  herein  runs  on  all  the  8  processors  of  MIT’s 
parallel  computer  Pleiades. 


5.1.4  Vertical  Flow  Velocity  Results 

The  model  is  run  with  grid  numbers  Nx  =  200,  Ny  =  200,  and  Nz  =  35,  i.e.,  simulat¬ 
ing  the  350-m  mixed  layer  observed  during  AUV  Mission  B9804107.  The  temperature 
and  salinity  profiles  during  AUV  Mission  B 9804 107  are  used,  as  shown  in  Figure  9.8. 
Model  computation  outputs  are  temperature  T,  horizontal  flow  velocities  u  and  v. 
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Table  5.2:  Physical  parameters 


Parameter 

Symbol 

Value 

Coriolis  frequency 

/ 

Coefficient  of  thermal  expansion 

a 

I 

Coefficient  of  saline  contraction 

P 

■rHRHKlSM 

Specific  heat  capacity 

cp 

Horizontal  diffusivity  of  momentum  and  temperature 

Vh,  Kh 

0.1  m2  s  1 

Vertical  diffusivity  of  momentum  and  temperature 

Uy  ,  Ky 

Acceleration  of  gravity 

9 

At  the  latitude  of  AUV  Mission  B9804107  site,  (p  fa  57°N,  so  we  have 
*  /  =  2 flsin((f))  fa  1.2  x  10-4  rad  s_1  where  ft  =  27r/day  is  the  angular 
velocity  of  Earth’s  rotation. 


and  vertical  flow  velocity  w. 

At  26520  s  (about  7.4  hours)  after  surface  cooling  starts,  the  vertical  flow  velocity 
w  is  shown  in  Figure  5.2.  The  upper  panel  displays  the  horizontal  cross-section 
at  the  250-m  depth,  the  same  depth  of  AUV  Mission  B9804107.  The  lower  panel 
displays  the  vertical  cross-section  at  y  =  1000  m.  Convective  cells  with  periodicity 
of  200  m  ~  250  m  are  observable  in  both  panels,  which  is  close  to  the  theoretical 
prediction  in  Equation  (5.3).  As  a  “cut”  out  of  the  lower  panel  of  Figure  5.2,  a  one¬ 
dimensional  plot  of  w  at  depth  250  m  is  shown  in  Figure  5.3.  Based  on  the  model 
output  for  two  hours  (from  5.4  hours  to  7.4  hours  after  the  onset  of  surface  cooling), 
we  compute  the  temporal-spatial  PSD  of  convective  vertical  velocity  w,  as  will  be 
presented  in  Subsection  6.2.1. 

5.2  Internal  Waves 

5.2.1  Background 

Internal  waves  occur  in  the  ocean’s  interior.  It  is  the  water’s  response  to  a  disturbance 
to  its  equilibrium  of  hydrostatically  stable  density  stratification,  via  the  gravitational 
restoring  force  [40],  [58],  [35].  Unlike  convection  which  occurs  in  a  vertically  mixed 
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W  (m/s).  Time  step  2652.  Initially  experimental  profiles. Depth  250m. 


-0.15  -0.1  -0.05  0  0.05 


Y=t000m. 


-0.15  -0.1  -0.05  0  0.05 


Figure  5.2:  Horizontal  and  vertical  cross-sections  of  vertical  flow  velocity  w  at  time 
7.4  hours. 

(i.e.,  unstratified)  water  column,  internal  waves  are  found  in  stably  stratified  water. 
Stable  stratification  is  depicted  by  the  buoyancy  frequency  N  (also  called  the  Brunt- 
Vaisala  frequency)  [9] ,  [59] : 


where  g  is  the  acceleration  of  gravity;  p(z )  is  the  water  density  as  a  function  of  z;  ,p0 
is  the  reference  density.  Note  that  z  points  upwards. 

Internal  waves  play  an  important  role  in  mass  and  momentum  transfer  in  the 
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W  (nVs).  Time  step  2652.  Y=  1000m.  Depth  250m. 


Figure  5.3:  w  v.s.  x  at  depth  250  m,  y=1000  m,  and  time  7.4  hours. 

ocean  [40].  Their  dynamics  is  essential  for  understanding  the  ocean  circulation  and 
the  temperature  and  salinity  structure  [40].  In  another  aspect,  the  sound-speed  fluc¬ 
tuations  induced  by  internal  waves  are  the  dominant  sources  of  the  high-frequency 
variability  of  acoustic  wave  fields  in  the  ocean  [60] . 

Vertical  flow  velocity  is  a  key  quantity  and  signature  of  ocean  convection  and  in¬ 
ternal  waves  [47],  [16].  For  the  purpose  of  experimental  test,  an  AUV-borne  Doppler 
sonar  acquired  flow  velocity  measurement  during  the  1998  Labrador  Sea  Experiment, 
as  will  be  presented  in  Chapter  9.  Thus  we  select  vertical  flow  velocity  as  the  quantity 
to  use  for  distinguishing  ocean  convection  from  internal  waves.  Prospects  of  intro¬ 
ducing  more  quantities  (e.g.,  temperature)  to  improve  classification  will  be  discussed 
in  Section  10.3. 

5.2.2  Vertical  Flow  Velocity  Spectrum  Based  on  Garrett- 
Munk-79  Model 

The  internal  wave  spectral  model  presented  by  C.  Garrett  and  W.  Munk  [2],  [61],  [40], 

[62] ,  the  so-called  GM79  model,  is  widely  used  in  the  oceanographic  community  [60], 

[63] ,  [64].  The  model  gives  the  power  spectrum  of  vertical  displacement  at  each  mode 
i  ( i  >  1): 


77 


Sc(uj,i) 


2b2 E0  N0  (i2  +  z2)  1  lo0(cj2  —  cog)1/2 
*  N{z)  E,“i[(i2  +  i2)"1]  ^ 


(5.17) 


where  N(z)  is  the  profile  of  buoyancy  frequency;  N0  is  the  surface  extrapolated  buoy¬ 
ancy  frequency;  b  is  the  e-folding  depth  of  N(z)\  uj0  is  the  Coriolis  frequency;  E0  is 
the  power  parameter;  z*  is  the  mode  scale  number.  The  values  of  these  quantities  in 
our  computations  will  be  given  in  Table  5.3. 

The  spectrum  of  vertical  velocity  W  can  be  deduced  from  Equation  (5.17)  by  the 
differentiation  relation 


(5.18) 


and  we  obtain 


Sw(w,i)  =  \joj\2S$(u,  i) 

=  u2S((u),  i) 

__  2 b2 E0  N0  (z2  +  z2)"1  a>0(u;2  —  Wg)1/2 

tt  N(Z)Y:z1[^+ii)-1}  z 


(5.19) 


which  follows  the  same  approach  as  in  [61]  whereby  the  horizontal  velocity’s  spectrum 
is  deduced  from  that  of  the  horizontal  displacement. 

In  the  thesis,  we  are  concerned  about  the  temporal-spatial  PSD  S'w  radial(u>,k) 
where  u  is  the  frequency  and  k  is  the  horizontal  wavenumber.  Note  that  based  on  the 
isotropy  assumption  [62],  only  a  single  radial  horizontal  wavenumber  k  =  y/k\  +  k 2 
is  used.  To  transform  Sw{u,i)  to  S'w_radial{u,  k)  (footnote  “radial”  is  added  for 
emphasis),  we  resort  to  the  dispersion  relation  [62]: 
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m  =  k 


(5.20) 


lN{zf 


.7 r  N{z)1 2-u 2 * 


UT  —  UJn 


N2 


or 


where  m  is  the  vertical  wavenumber.  It  leads  to  the  k  —  i  relation 


to* 


(5.21) 


which  enables  the  transformation 


-radial  0*5  ^0  (^5 


(5.22) 


which  is  equivalent  to 


-radial  0*5  ^ 


.  7T  U2  —  U>? 


b  V  N2  -  U 


7i)  =  (z\h^j)Sw(u,i)  (5.23) 


7T  v  w*  -  u;£ 


Evaluation  of  Equation  (5.23)  is  carried  out  numerically.  The  procedure  is  as 
follows: 


1.  At  any  specified  u,  map  the  wavenumber  range  K  =  [kmin,kmax]  to  the  mode 
number  range  I  using  Equation  (5.21).  Note  that  the  mode  number  must  be 
an  integer.  So  in  the  mapping,  the  mode  number  range  /  is  obtained  after 
rounding.  The  spectrum  is  then  computed  directly  by  Equation  (5.19).  As 
shown  in  Equation  (5.23),  the  scaling  ensures  preservation  of  power  in  different 
spaces. 

2.  Due  to  the  rounding  in  the  above  step,  the  mapped  mode  number  range  I  does 

not  exactly  correspond  to  the  required  wavenumber  range  K,  but  to  a  varied 

range  K' .  Spectrum  values  on  K  is  obtained  by  interpolation  using  values  on 
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K'.  Note  that  Equation  (5.21)  dictates  that  at  any  specified  w,  k  has  a  lower 
bound  corresponding  to  mode  number  2=1: 


Slower  -bound  — 

below  which  S'W  radial(u,  k)  is  assigned  zero. 

3.  Loop  through  all  u  values  between  the  Coriolis  frequency  u>o  and  the  buoyancy 
frequency  N(z). 

Hereafter  we  use  frequency  and  wavenumber  notations  77  =  and  uradia,  =  Jl. 

ZTT  2l T 

Correspondingly,  the  PSDs  using  the  two  sets  of  notations  are  expressed  as 


I UJ2 

\In! 


jradial^i^)  ‘“’W -radiali^^Vi  ‘^TTVradial')  &W -radial (Vi  ^radial) 


(5.25) 


Table  5.3:  GM79  model  parameters 


N0  (rad/s) 

N(z )  (rad/s) 

b  (m) 

Eq 

Imax 

cao  (rad/s) 

5.2  x  10“3 

5.0  x  10~3 

1.3  x  103 

6.3  x  10“5 

3 

100 

1.2  x  10-4 

In  the  computation  of  Equation  (5.19)  and  Equation  (5.23),  the  following  pa¬ 
rameter  values  are  taken  from  [62]:  5=1300  m,  E0  =  6.3  x  10“5  (dimensionless), 
N0  =  5.2  x  10  3  rad/s  (equivalent  to  about  three  cycles  per  hour),  i*=3.  Those 
parameters  have  also  been  adopted  in  other  research  [60].  The  two  remaining  pa¬ 
rameters,  the  Coriolis  frequency  loq  and  the  buoyancy  frequency  7V(z)  are  assigned 
in  the  thesis  context.  In  relation  to  the  Labrador  Sea  Experiment  data  (to  be  pre¬ 
sented  in  Chapter  9),  we  use  latitude  4>  =  5T°N  of  the  experimental  site  in  calculating 
coq  =  2fl  sinif)  «  1.2  x  10  4  rad/s  where  =  27r/day  is  the  angular  velocity  of  Earth’s 
rotation.  We  let  N(z)  =  5.0  x  10-3  rad/s  be  very  close  to  the  upper  bound  jV0. 
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In  Chapter  7,  we  will  look  at  the  internal  wave  spectrum  corresponding  to  a  lower 

N(z). 

The  calculated  Sw^adiaiiv,  "radial)  is  shown  in  Figure  5.4.  As  expressed  by  Equa¬ 
tion  (5.24),  at  any  cj,  k  has  a  lower  bound  corresponding  to  mode  number  i  =  1,  below 
which  S'w_radial(u,  k)  vanishes.  This  is  why  we  see  a  “shadow  region”  in  Figure  5.4, 
and  equivalently  a  “shadow  disk”  in  Figure  5.5. 


10l°9iOSw.radia|(7l’VHm3s'1) 


Figure  5.4:  Internal  wave  vertical  velocity’s  temporal-spatial  (radial)  PSE) 

Sw-radial{j]i  ^radial)-  '  - 

"  Using  Equation  (A. 5)  and  Equation  (A.8)in  Appendix  A,  the  “for  line  survey” 
PSD  Sw(v ,  ^)  can  be  obtained  by  an  integration  over  one  of  the  two  wavenumbers, 
as  illustrated  in  Figure  5.5.  Sw{n ,v)  is  shown  in  Figure  5.6  (only  its  first  quad¬ 
rant  is  displayed).  Note  that  after  integration  over  one  wavenumber  (as  expressed 
by  Equation  (A.5)),  the  shadow. region  as  seen  in  Figure  5.4  no  longer  exists.  This 
phenomenon  can  be  explained  by  the  illustration  of  Figure  5.5. 

By  integrating  Sw{v>'v)  over  z/  or  77,  we  can  obtain;  the  vertical  velocity’s  77- 
spectrum  and  ^-spectrum,  as  shown  in  the  upper  and  lower  panels  of  Figure  5.7, 


kVl 


Isotropic  contours  of 

S3d(T|,Vi,V2) 


Integration  line  for 
producing  S(r|,vi) 


"Shadow”  region 

at  r|=co/(27i) 


Vi 


Figure  5.5:  Integration  over  one  wavenumber  to  produce  the  “for  line  survey”  PSD 
Sw(v,  u)  °f  internal  wave  vertical  velocity. 


respectively.  For  the  sake  of  comparison,  those  two  spectra  are  also  calculated  for 
the  vertical  displacement,  as  shown  in  Figure  5.8.  Due  to  the  differentiation  relation 
as  expressed  in  Equation  (5.18),  vertical  velocity’s  r/-spectrum  is  flattened  up  from 
vertical  displacement’s  r]~2  asymptote. 
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1Oi°910Sw(ll'v)  (m3s  ^ 


Figure  5.6:  Internal  wave  vertical  velocity’s  temporal-spatial  (for  line  survey)  .  PSD 
SwiViV)  (the  first  quadrant). 


Frequency  power  spectrum  density  of  vertical  velocity 


Figure  5.7:  Internal  wave  vertical  velocity’s  frequency-spectrum  (upper)  and 
wavenumber-spectrum  (lower). 
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Wavenumber  power  spectrum  density  of  vertical  displacement 


Figure  5.8:  Internal  wave  vertical  displacement’s  frequency-spectrum  (upper)  and 
wavenumber-spectrum  (lower). 


Chapter  6 


AUV-Based  Spectral  Classification 
of  Ocean  Convection  and  Internal 
Waves 


6.1  Classification  Steps 

We  now  apply  the  AUV-based  classifier  to  distinguish  ocean  convection  from  internal 
waves.  According  to  the  classifier  architecture  illustrated  in  Figure  4.1,  the  steps 
of  classification  are  summarized  as  follows  (internal  wave  is  denoted  class  1  while 
convection  is  denoted  class  2): 

1.  For  each  process,  derive  the  mingled  spectrum  SY(f )  from  its  temporal-spatial 
spectrum  Sxijhv)  by  Equation  (3.4).  Then  by  including  the  data  window’s 
effect,  calculate  class  mean  vectors  MYl  and  MY2  by  Equation  (4.20)  and  Equa¬ 
tion  (4.21),  respectively. 

2.  Obtain  the  transformation  vector  V  by  Equation  (4.25). 

3.  By  Equation  (4.26),  we  apply  V  to  the  PSD  estimate  SY(f)  of  the  input  time 
series  Y(t).  The  resultant  scalar  feature  z  is  compared  with  a  threshold  to  make 
the  classification  decision. 
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In  this  chapter,  we  use  model-based  synthesized  data  in  step  3.  We  evaluate  the 
classifier’s  performance  by  the  statistics  of  the  ensemble  of  z.  In  Chapter  9,  we  will 
use  experimental  data  to  test  the  classifier. 

6.2  Mingled  Power  Spectrum  Density  (PSD)  of 
Vertical  Velocity  of  Convection  and  Internal 
Waves 

6.2.1  Temporal-Spatial  PSDs  of  the  Two  Processes 

1.  Ocean  convection. 

In  Section  5.1.4,  we  have  configured  and  run  a  2000  m  x  2000  m  x  350  m  mixed- 
layer  convection  model,  over  a  time  duration  of  two  hours.  The  model  has  a  grid 
size  of  10  m  and  a  data  output  step  of  30  s.  At  any  depth,  there  are  thus  200 
lines  evolving  for  240  time  steps.  We  carry  out  model  computations  at  the  depth 
of  250  m,  selected  to  coincide  with  that  of  the  AUV’s  Labrador  Sea  experimental 
data  which  will  be  presented  in  Section  9.4  and  tested  in  Section  9.5. 

To  obtain  vertical  velocity’s  temporal-spatial  PSD  for  a  line  AUV  survey,  we 
conduct  two-dimensional  Fourier  transform  [65]  (over  time  and  distance)  for 
each  line.  Then  we  take  the  average  of  the  200  lines’  two-dimensional  peri- 
odograms  as  the  averaged  periodogram.  This  procedure  is  illustrated  by  Fig¬ 
ure  6.1.  While  the  mingled  spectrum  principle  does  not  rely  on  isotropy,  we 
consider  the  convection  field  to  be  isotropic  because  of  the  isotropic  surface 
cooling  in  the  model.  Based  on  symmetry  properties  as  expressed  by  Equa¬ 
tion  (A.  10)  and  Equation  (A.  12),  the  whole  tj-v  spectrum  is  constructed  using 
the  first  quadrant  of  the  averaged  periodogram.  The  resultant  spectrum  is 
shown  in  the  upper  panel  of  Figure  6.2.  It  is  regarded  as  the  temporal-spatial 
PSD  template  for  convection. 
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Temporally,  the  vertical  velocity  field  varies  little  during  the  two  hours  of  evolu¬ 
tion  (from  5.4  hours  to  7.4  hours  after  the  onset  of  surface  cooling)  as  convection 
approaches  a  stationary  state.  The  observed  baseband  spectrum  on  the  77-axis 
is  mostly  due  to  the  two-hour  window  effect,  on  the  zz-axis,  however,  there  is  a 
peak  at  about  0.005  m-1  because  convective  cells  have  a  spatial  periodicity  of 
about  200  m  (as  displayed  in  Figure  5.2).  As  will  be  shown  in  Subsection  6.2.2, 
we  can  utilize  the  AUV’s  speed  to  highlight  this  feature  of  convection  for  clas¬ 
sification  against  internal  waves. 


Sx(t|,v)  template. 

Figure  6.1:  Generation  of  the  temporal-spatial  PSD  of  convective  vertical  velocity 
from  model  data. 

2.  Internal  waves. 

The  temporal-spatial  PSD  of  internal  wave’s  vertical  velocity  has  been  obtained 
in  Subsection  5.2.2,  as  shown  in  Figure  5.6.  That  spectrum  is  confined  within  an 
upper  bound  of  buoyancy  frequency  {N(z)  =  8  x  10“4  Hz,  equivalent  to  nearly 
three  cycles  per  hour)  on  the  77-axis.  In  the  ocean,  however,  processes  with 
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Figure  6.2:  Temporal-spatial  PSD  of  vertical  velocity  of  convection  (upper)  and  in¬ 
ternal  waves  (lower,  with  extended  plateau). 

frequencies  higher  than  the  buoyancy  frequency  do  exist  [40],  like  turbulence  [9]’. 
We  therefore  need  to  consider  higher-frequency  processes  along  with  internal 
waves. 

According  to  the  Power  Spectrum  Density  (PSD)  of  ocean  wave  kinetic  en- 
•  ergy  as  illustrated  in  Figure  6.3,  the  spectrum  follows  r]~2  within  the  internal 
wave  frequency  range  ,  (i.e.,  between  the  Coriolis  frequency  and  the  buoyancy 
frequency)  and  follows  ry”4  above  the  buoyancy  frequency..  The  observed  rj~2 
power  law  is  consistent  with  observations  in  classic  internal  wave  papers  [2],  [61], 


Figure  6.3:  Notional  plot  of  frequency  spectrum  of  kinetic  energy  of  internal  wave 
and  higher-frequency  waves.  The  turning  point  is  at  the  buoyancy  frequency.  (Based 
on  a  plot  with  courtesy  of  Dr.  Thomas  Curtin.) 


[40],  [62],  In  Figure  6.3,  the  frequency  bounds  are:  Coriolis  frequency=10~5  Hz, 
corresponding  to  a  latitude  of  about  30o1;  buoyancy  frequency=2.8  x  10~4  Hz, 
equivalent  to  about  one  cycle  per  hour;  upper  end  of  frequency=0.017  Hz,  up  to 
where  the  r/“4  power  law  still  applies  well.  By  integration,  the  power  contained 
above  the  buoyancy  frequency  is  found  to  be  about  1%  of  the  internal  wave 
power.  Varying  the  buoyancy  frequency  up  to  three  cycles  per  hour  (as  used  for 
our  internal  wave  spectrum  computation  so  far),  this  ratio  rises  to  about  3%. 
To  account  for  power  contribution  from  higher-frequency  processes,  we  need  to 
extend  internal  wave’s  power  spectrum  to  above  the  buoyancy  frequency.  Before 
carrying  out  the  extension,  we  should  determine  the  upper  bounds  of  frequency 


1  At  a  higher  latitude  (e.g.,  ~  57°N  for  the  1998  Labrador  Sea  Experiment  to  be  presented 

in  Chapter  9),  the  Coriolis  frequency  will  be  higher.  The  ratio  above  huowincy  frequency 

...  ,.  ,  Aii.-  ,  .  .  .  .  ,  Power  of  internal  waves  , 

will  correspondingly  rise.  As  to  be  given  below,  the  ratio  is  set  large  enough  to  allow  for  this 

variation. 
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and  wavenumber: 


•  Temporal  frequency  rj 

In  convective  vertical  velocity  PSD  computation,  the  30-s  time  step  in 
running  the  model  determines  that  the  upper  end  of  temporal  frequency 
is  I  x  30s  =  0-0167  Hz.  For  internal  wave  vertical  velocity,  most  power 
is  contained  below  0.0167  Hz  which  corresponds  to  the  highest  frequency 
point  drawn  in  Figure  6.3.  So  we  take  0.0167  Hz  as  the  upper-end  temporal 
frequency  for  both  convection  and  internal  wave  PSDs. 

•  Spatial  frequency  v 

In  convective  vertical  velocity  PSD  computation,  the  10-m  grid  size  in 
running  the  model  determines  that  the  upper  end  of  spatial  frequency  is 
2  x  10m  =  6-05  m-1.  For  internal  wave  vertical  velocity  PSD  as  shown 
in  Figure  6.3,  the  power  density  is  very  low  above  0.05  m-1  (power  above 
this  spatial  frequency  accounts  for  less  than  0.3%  of  the  total  power).  So 
we  take  0.05  m_1  as  the  upper-end  spatial  frequency  for  both  convection 
and  internal  wave  PSDs. 

Now  for  the  temporal-spatial  PSD  of  internal  wave  vertical  velocity,  we  add  a 
spectrum  plateau  above  the  buoyancy  frequency  to  account  for  higher-frequency 
processes.  Although  a  plateau  is  not  an  accurate  description  of  the  spectrum, 
we  deem  it  sufficing  to  serve  the  purpose  of  thesis  work  since  the  forthcoming 
computation  of  mingled  spectrum  is  in  an  integration  sense.  From  the  per¬ 
spective  of  classification,  the  spectrum  extension  will  prevent  a  classifier  from 
unduly  taking  advantage  of  a  vanishing  part  of  a  spectrum2. 

Calculated  based  on  Figure  6.3,  the  ratio  frequency  ,s 

°  ’  Power  ot  internal  waves 

used  to  set  the  height  of  plateau.  Although  the  ratio  calculated  this  way  is  for 
kinetic  energy  but  not  for  (vertical  velocity)2  per  se,  we  still  deem  it  a  usable 

2 In  detection  theories,  the  problem  is  referred  to  as  “singular  detection”  [21],  [66]. 
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reference.  To  allow  for  unaccounted  power  at  even  higher  frequency  and  also 
allow  for  latitude  variations,  we  set  the  power  ratio  to  0.1. 

With  this  plateau  extension  above  the  buoyancy  frequency,  the  77  —  v  spectrum 
for  internal  wave  vertical  velocity  is  shown  in  the  lower  panel  of  Figure  6.2. 
Symmetries  are  due  to  isotropy,  as  expressed  by  Equation  (A. 10)  and  Equa¬ 
tion  (A. 12).  Note  that  for  the  sake  of  testing  the  classifier,  we  have  scaled 
internal  wave  vertical  velocity’s  amplitude  such  that  its  power  equals  that  of 
convective  vertical  velocity.  On  the  77-axis,  internal  wave  power  is  confined  be¬ 
tween  the  Coriolis  frequency  (1.2  x  10“4  rad/s,  corresponding  to  the  Labrador 
Sea  Experiment  latitude)  and  the  buoyancy  frequency  (5.0  x  10~3  rad/s,  i.e., 
nearly  three  cycles  per  hour).  On  the  z/-axis,  most  power  lies  at  very  low 
wavenumber,  showing  no  peak  at  away  from  v  -  0.  This  distinction  from 
convective  vertical  velocity  is  what  a  cruising  AUV  can  take  advantage  of  for 
classification. 

6.2.2  Mingled  PSDs  at  a  Series  of  AUV  Speeds 

Having  obtained  Sx{rj,u)  of  convective  and  internal  wave  vertical  velocities,  let  us 
derive  the  AUV-seen  SY{f)  by  the  mingled  spectrum  principle,  as  illustrated  in  Fig¬ 
ure  6.4.  We  intend  to  depict  the  temporal-spatial  distinctions  of  the  two  processes 
in  Figure  6.4,  rather  than  accurately  plot  Sx{t],v).  Convection’s  spatial  peak  on 
the  I'-axis  is  projected  onto  the  /-axis  of  the  corresponding  mingled  spectrum.  At  a 
higher  AUV  speed,  the  spectral  peak  on  the  /-axis  is  pulled  farther  away  from  /  =  0. 
For  internal  waves,  however,  the  picture  is  different.  Because  the  internal  wave  power 
is  concentrated  at  baseband  on  the  77-axis  and  the  ^-axis,  the  corresponding  mingled 
spectrum  also  lies  at  baseband  on  the  /-axis.  A  higher  AUV  speed  will  not  change 
this  basic  spectral  shape.  Based  on  this  inspection  even  before  conducting  computa¬ 
tions,  we  project  that  the  classifier’s  performance  will  improve  with  AUV  speed  as 
convection’s  spatial  feature  is  highlighted  in  contrast  with  internal  waves. 

We  apply  Equation  (3.4)  at  a  series  of  AUV  speeds  u=  1  m/s,  0.25  m/s,  0.1  m/s, 


92 


Internal  Wave  Convection 


Figure  6.4:  Derivation  of  mingled  spectra  from  rj-v  spectra  of  convection  and  internal 
wave  vertical  velocities. 

and  0.05  m/s.  The  resultant  mingled  spectra  of  convective  and  internal  wave  vertical 
velocities  are  compared  in  Figure  6.5.  The  effect  of  data  window  has  been  included  in 
the  calculations,  by  using  Equation  (4.20)  and  Equation  (4.21).  So  the  results  shown 
in  Figure  6.5  are  class  mean  vectors  Myi_o(f)  and  My2_o(/)-  The  data  window 
length  is  set  to  1400  s  to  coincide  with  that  of  the  AUV’s  Labrador  Sea  experimental 
data  which  will  be  presented  in  Section  9.4  and  tested  in  Section  9.5.  The  results 
in  Figure  6.5  are  consistent  with  the  predictions  inspected  from  Figure  6.4:  a  higher 
AUV  speed  improves  classification  by  highlighting  convection’s  spatial  feature  through 
mingled  spectrum  projection. 
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AUV  speed  =  1  m/s,  with  1400-s  window  effect. 


-0.02  -0.015  -0.01  -0.005  0  0.005  0.01  0.015  0.02 

Figure  6.5:  Mingled  spectra  of  vertical  velocity  of  convection  and  internal  waves!  ' 


6.3  Feature  Transformation  Vector 

Having  obtained  class  mean  vectors  Myi_0(/)  and  MY2_q{f),  the  feature  transforma¬ 
tion  vector  V  is  derived  by  Equation  (4.25).  For  an  AUV  speed  of  1  m/s,  V  is  shown 
in  the  lower  panel  of  Figure  6.6,  labeled  “unmodified”.  In  Chapter  7,  the  computation 
method  of  V  will  be  modified  in  consideration  of  parameter  uncertainties.  Although 
we  do  not  present  the  modification  method  herein  for  the  sake  of  clarity,  we  show  the 
modified  V  and  the  resultant  classifier  performance  hereafter.  Class  mean  vectors 
Afyi_o(/)  and  My2_0(/)  and  the  transformation  vector  V  for  AUV  speed  series  u— 
1  m/s,  0.25  m/s,  0.1  m/s,  and  0.05  m/s,  are  shown  in  Figure  6.7. 

For  the  purpose  of  the  forthcoming  experimental  data  test,  we  have  normalized  the 
power  of  Myi  o  and  My2_o  to  that  of  the  experimental  data  set  which  will  be  presented 
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AUV  speed  =  1  m/s 


Figure  6.6:  Class  mean  vectors  Myi_0  and  My2_ o,  and  the  feature  transformation 
vector  V  for  AUV  speed  1  m/s. 

in  Section  9.4.  The  upper  frequency  bound  is  prescribed  such  that  both  Myi_o  and 
My2J)  are  higher  than  0.0015  (m/s)2/Hz.  This  power  density  value  maintains  an  SNR 
of  about  20  dB  over  the  instrument  noise  floor  of  the  Acoustic  Doppler  Velocimeter 
(ADV)  under  normal  operation  conditions  (an  AUV-borne  ADV  acquired  the  flow 
velocity  data  in  the  Labrador  Sea  Experiment,  as  will  be  presented  in  Chapter  9).  At 
a  lower  AUV  speed,  levels  of  mean  spectra  decrease  at  high  frequency,  thus  the  valid 
frequency  range  shrinks  as  the  vehicle  speed  drops. 

In  computations  presented  in  this  chapter,  we  let  vejf  =  1  (please  refer  to  Equa¬ 
tion  (4.25)).  This  implies  that  no  time-domain  segmentation  or  frequency-domain 
smoothing  is  applied  in  periodogram  computation.  The  periodogram  thus  generated 
is  noisy.  For  the  AUV-based  classification,  treating  this  adverse  situation  is  necessi¬ 
tated  when  the  data  length  is  small,  since  segmentation  would  lead  to  an  even  shorter 
duration  for  each  segment,  which  corresponds  to  a  wider  bandwidth  Bw.  In  the  formu¬ 
lation  of  the  classifier’s  transformation  vector  V,  it  is  required  that  adjacent  frequency 
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Figure  6.7:  Class  mean  vectors  MY i_0  and  My2_o  and  the  feature  transformation 
vector  V  for  a  series  of  AUV  speeds. 


points  are  separated  by  a  spacing  of  Bw  such  that  the  PSD  estimates  are  uncorrelated. 
Given  a  limited  frequency  range  constrained  by  a  signal-to-noise-ratio  (SNR)  over 
the  instrument  noise  floor,  the  total  number  of  valid  frequency  points  participating 
in  classification  will  decrease  when  Bw  increases.  Frequency-domain  smoothing,  on 
the  other  hand,  also  reduces  the  total  number  of  valid  frequency  points  participating 
in  classification,  because  smoothing  makes  the  PSD  estimates  within  the  smoothing 
band  no  longer  uncorrelated.  The  Labrador  Sea  AUV  data  used  in  the  thesis  has  a 
length  of  about  1400  s,  corresponding  to  Bw  «  7  x  10-4  Hz.  With  no  time-domain 
segmentation  or  frequency-domain  smoothing  applied  in  periodogram  computation, 


i.e,  ueff  =  1,  the  total  number  of  valid  frequency  points  participating  in  classification 
are:  14,  9,  6,  5  for  AUV  speed  1  m/s,  0.25  m/s,  0.1  m/s,  and  0.05  m/s,  respectively,  as 
shown  in  Figure  6.7.  Hence  further  segmentation  or  smoothing  would  be  unrealistic 
for  this  case. 

An  AUV  should  be  equipped  with  the  ability  of  carrying  out  classification  based 
on  a  short  data  record.  This  is  necessary  for  the  vehicle  to  make  timely  responses  to 
detected  processes  of  interest.  On  the  other  hand,  the  classifier’s  underlying  assump¬ 
tions  of  temporal  stationarity  and  spatial  homogeneity  would  become  more  challenged 
when  the  data  length  is  large.  With  the  above  considerations,  we  deem  it  necessary 
to  investigate  the  classifier’s  performance  under  veff  =  1  which  constitutes  the  most 
challenging  situation  in  the  respect  of  available  data  length.  This  is  also  necessary  for 
testing  the  classifier  with  the  short  experimental  data.  In  Chapter  7,  we  will  present 
the  results  for  veff  =  4  by  assuming  a  longer  data  record. 

6.4  Model-Based  Simulations 

6.4.1  Simulation  of  AUV-Acquired  Convective  Vertical  Ve¬ 
locity 

For  simulating  a  line  survey  in  the  convection  field,  the  AUV-recorded  time  series 
are  directly  drawn  from  the  convection  model  at  depth  250  m  which  is  displayed 
in  Figure  5.2.  We  thus  consider  each  of  the  200  lines  to  be  an  AUV  survey  line  with 
a  duration  of  no  more  than  two  hours,  as  illustrated  in  Figure  6.1.  At  an  AUV  speed 
<  |  m/s,  interpolation  is  done  in  distance  to  fit  the  initial  time  step  of  30  s;  at  an 
AUV  speed  >  |  m/s,  interpolation  is  done  in  time  to  fit  the  initial  grid  size  of  10  m. 
An  example  of  an  AUV-recorded  vertical  velocity  time  series  is  shown  in  Figure  6.8. 
In  the  second  panel,  circles  mark  the  original  data  points,  and  crosses  mark  the 
interpolated  data. 

An  ensemble  of  200  AUV  survey  lines  are  used  for  a  classifier  test  at  each  prescribed 
AUV  speed.  To  add  randomness  to  different  test  runs,  the  starting  time  and  location 
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kline=45,  time  starts  at  60  s  and  location  starts  at  240  m. 


Vertical  velocity  time  series 


Figure  6.8:  An  AUV  survey  line  in  the  convection  field  at  speed  1  m/s. 

are  randomly  picked  (within  range).  For  example,  in  the  test  run  corresponding 
to  Figure  6.8,  the  start  time  is  60  s  while  the  start  location  is  240  m. 

6.4.2  Simulation  of  AUV-Acquired  Internal  Wave  Vertical 
Velocity 

At  AUV  speed  1  m/s,  the  mingled  PSD  of  internal  wave  vertical  velocity  SY_intwave(f ) 
(obtained  by  Equation  (3.4))  is  shown  by  the  solid  curve  in  the  second  panel  of  Fig¬ 
ure  6.9  (without  the  data  window’s  effect).  SY_jntwave(f)  approximately  follows  a 
power  law  of  f°  at  low  frequency,  and  follows  /~2  at  high  frequency.  This  property 
reminds  us  of  the  feasibility  of  simulating  the  corresponding  time  series  as  the  output 
of  a  first-order  AutoRegressive  ( AR)  model  [67] . 

The  difference  equation  for  a  first-order  AR  model  is 
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Mingled  PSD  for  internal  wave  vertical  velocity  at  AUV  speed  1  m/s 


Figure  6.9:  The  mingled  PSD  of  internal  wave  vertical  velocity  v.s.  the  AR  modeled 
PSD  at  AUV  speed  1  m/s. 


y(n )  =  ay(n  -  1)  +  u(n)  (6.1) 

where  a  is  an  coefficient  constrained  by  |a|  <  1;  u{n)  is  a  white  noise  process  with 
variance  o\. 


rc 

output  y(n) 

white  noise  u(n)  V 

r 

— « —  z'1  — - 

Figure  6.10:  A  first-order  AR  model. 


Consider  u(n )  and  y(n)  to  be  the  input  and  output  of  a  system,  and  then  the 
system’s  transfer  function  is 


99 


(6.2) 


H(z) 


Y(z) 

U(z) 


1 

1  —  OLZ~X 


which  is  depicted  by  Figure  6.10.  The  coefficient  a  represents  the  pole  z0  =  a.  Its 
closeness  to  the  unit  circle  determines  the  frequency  response  of  the  system  H(z). 

Since  the  input  u(n)  is  white  noise  of  power  o\,  the  PSD  of  the  output  y(n)  is 
simply 


SY-AR(f)  =  al  \H{f)\2  = 


at 


cr: 


|1  —  ae  |2  1  —  2a  cos(2nf )  +  a2 


(6.3) 


To  facilitate  a  quick  observation  of  SY_AR(f ),  we  can  utilize  the  technique  of  Bode 
diagram.  The  denominator  of  Sy_ar(I )  in  Equation  (6.3)  can  be  written  as 


|1  —  ae  j27r^|2  =  [1  —  a  cos(2nf  )]2  +  [ asin(2irf )]2  (6.4) 

asin(- A-)—? 

At  /o  =  - - - ,  1  —  acos(2ivf )  equals  asinfaf)  (condition:  a  >  ) .  For 

f  <  fo,  the  former  dominates;  for  f  >  f0,  the  latter  dominates.  As  a  Bode  diagram 
approximation,  we  only  keep  the  dominant  terms  below  or  above  /0.  For  /  c  1,  we 
can  further  apply  the  Taylor  series  expansions  to  sin(-)  and  cos(-)  functions  and  only 
keep  the  first  terms,  which  leads  to  f0  «  Equation  (6.3)  is  then  approximated 
by 


Sr(f) 


<  (1  — a)2  0  <  /  <  fo 

k  (4  fo  <  f  «  i 


(6.5) 


where  Sy(f)  follows  f°  below  /o  and  follows  f~ 2  above  /o,  respectively.  This  verifies 
the  applicability  of  simulating  the  AUV-acquired  time  series  of  internal  wave  vertical 
velocity  by  use  of  a  first-order  AR  model.  Figure  6.11  displays  the  comparison  of 


100 


an  AR  model  output  PSD  and  its  Bode  diagram  approximation.  The  first-order  AR 
model  design  procedure  can  be  summarized  as  follows: 


Figure  6.11:  Comparison  of  a  first-order  AR  modeled  PSD  and  its  Bode  diagram. 

1.  Inspect  the  Bode  diagram’s  /0  for  the  given  PSD.  Determine  an  initial  AR 
coefficient  a  by  use  of  fo  —  Assign  a  corresponding  value  to  white  noise 
variance  cr^  to  guarantee  that  the  AR  model  output  has  the  same  power  as  that 
of  the  given  PSD. 

2.  Do  a  finer  search  of  a,  until  the  AR  output  PSD  best  matches  the  given  PSD. 
The  metric  is  ||5y(/)  —  Sy^ji(f)\\2  (weighting  can  be  applied  based  on  the 
feature  transformation  vector  V).  Each  time  a  is  adjusted,  so  is  to  preserve 
power. 

With  parameters  a  =  0.9627  and  o\  =  1.1041  x  10~4,  comparison  of  the  AR 
modeled  PSD  and  the  objective  PSD  of  internal  wave  vertical  velocity  is  shown  in  Fig¬ 
ure  6.9.  Including  the  effect  of  the  1400-s  data  window,  Figure  6.9  updates  to  Fig¬ 
ure  6.12.  The  two  PSDs  appear  to  match  well.  One  sample  time  series  generated  by 
this  first-order  AR  model  and  its  PSD  estimate  are  shown  in  Figure  6.13. 

At  a  lower  AUV  speed  such  as  0.1  m/s,  the  mingled  spectrum  of  internal  wave 
vertical  velocity  shows  a  spectral  peak  at  non-zero  frequency.  The  non-zero  spectral 
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Figure  6.12:  With  a  1400-s  data  window,  the  mingled  PSD  of  internal  wave  vertical 
velocity  v.s.  the  AR  modeled  PSD  for  AUV  speed  1  m/s. 

peak  prompts  us  to  use  a  second-order  AR  model,  which  has  a  pair  of  conjugate  poles, 
thus  providing  an  additional  degree  of  freedom  for  the  location  of  peak  frequency  [67]. 
The  difference  equation  for  a  second-order  AR  model  is 


y(n)  =  oiXy(n  -  1)  +  a2y{n  -  2)  +  u(n)  (6.6) 

where  ax  =  2 r  cos(2w  fpeak)  (r  is  the  radius  of  the  pair  of  conjugate  poles  and  fpeak  is 
the  peak  frequency)  and  a2  =  —r2\  u(n )  is  a  white  noise  process  with  variance  a2. 
The  system’s  transfer  function  is 


H(z) 


Y(z) 

U(z) 


_ 1 _ 

1  —  ol\z~1  —  a2z~2 


(6.7) 
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A  sample  time  series  of  AR  model  output  at  AUV  speed  1  m/s 


Figure  6.13:  A  sample  time  series  and  its  PSD  generated  by  a  first-order  AR  model 
for  AUV  speed  1  m/s. 


which  is  depicted  by  Figure  6.14.  The  PSD  of  the  output  y(n )  is  [67] 


SY_AR(f)=al\H(f)\2 

_ _ ol _ 

_  |1  -  a\e-***t  -  a2e~j4*f\2 

= _ _ 

—  y'g  j2rr(/  /peafc)|2  j  2  j’ZnU+fpeak)  |2 

For  AUV  speed  0.1  m/s,  a  second-order  AR  model  is  designed  to  match  the 
mingled  PSD.  The  resultant  model  parameters  are:  ol\  =  1.9217,  ol2  =  —0.9235,  and 
o\  =  4.5521  x  10-7.  With  the  effect  of  a  1400-s  data  window,  the  modeled  PSD  is 
compared  with  the  mingled  PSD  in  Figure  6.15. 

One  sample  time  series  generated  by  this  second-order  AR  model  and  its  PSD 
estimate  are  shown  in  Figure  6.16.  Like  in  the  convection  field,  200  AUV  survey  lines 
in  the  internal  wave  field  are  randomly  generated  in  each  classifier  test. 
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Figure  6.14:  A  second-order  AR  model. 


6.5  Classifier  Test  Results 

Now  we  test  the  classifier  by  convection  and  internal  wave  data  generated  in  the 
preceding  section.  For  each  time  series  of  AUV  data,  its  PSD  estimate  is  converted 
to  a  scalar  feature  z  by  the  transformation  vector  V  following  Equation  (4.1).  200 
lines  of  AUV  data  in  the  convection  field  and  another  200  lines  in  the  internal  wave 
field  are  used  for  each  test,  as  summarized  in  Table  6.1.  The  classifier’s  performance 
is  evaluated  by  the  statistics  of  the  resultant  ensemble  of  2. 


Table  6.1:  Parameters  and  method  for  classifier  test  of  convection  versus  internal 
waves 


Test  parameters 

Simulation  method 

Internal  wave  vertical  velocity 

N(z)  =  5.0  x  10-3  rad/s 

No  =  5.2  x  10~3  rad/s 

200  time  series  by 

AR  modeling 

Convective  vertical  velocity 

surface  heat  flux=300  W /m2 
mixed  layer  depth=350  m 

200  time  series  by 

extracting  AUV  survey  lines  from 

the  convective  box  at  depth  250  m 

Corresponding  to  Figure  6.5,  we  test  the  classifier  at  a  series  of  AUV  speeds  u= 
1  m/s,  0.25  m/s,  0.1  m/s,  and  0.05  m/s.  The  AUV’s  Labrador  Sea  experimental  data 
was  acquired  at  a  speed  of  about  1  m/s.  We  therefore  present  the  case  for  AUV  speed 
1  m/s  with  more  details  than  the  other  speeds,  for  clarity  and  also  to  serve  Section  9.5. 
Corresponding  to  200  AUV  survey  lines  in  the  internal  wave  field  and  200  lines  in 
the  convection  field,  we  get  the  histograms  of  z  as  shown  in  Figure  6.17.  In  this 
test  run,  the  distributions  of  z  in  the  two  classes  do  not  overlap.  The  corresponding 
Receiver  Operating  Characteristic  (ROC)  curve  is  shown  in  Figure  6.18  (please  refer 
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Mingled  PSD  (with  1400-s  window  effect)  for  internal  wave  vertical  velocity  at  AUV  speed  0.1  m/s 


Figure  6.15:  With  a  1400-s  data  window,  the  mingled  PSD  of  internal  wave  vertical 
velocity  v.s.  the  AR  modeled  PSD  for  AUV  speed  0.1  m/s. 

back  to  Section  4.4  for  explanations  of  ROC). 

The  arrow  in  Figure  6.17  marks  the  z  value  corresponding  to  the  AUV’s  Labrador 
Sea  data.  So  we  see  that  the  experimental  z  falls  in  the  cluster  of  the  model-based 
simulations.  We  will  present  the  experimental  data  in  Section  9.4  and  discuss  their 
tests  on  the  classifier  in  Section  9.5. 

With  the  same  procedure,  we  obtain  the  classifier’s  performance  at  other  AUV 
speeds.  For  the  speed  series,  Figure  6.18  shows  the  histograms  and  the  correspond¬ 
ing  ROC  curves.  The  corresponding  first  and  second-order  statistics  of  z  is  shown 
in  Table  6.2. 

Classification  of  the  two  processes  is  very  difficult  at  low  AUV  speeds  0.05  m/s  and 
0.1  m/s,  shown  by  the  severely  overlapping  histograms  and  the  corresponding  ROC 
curves  in  Figure  6.18.  As  noted  in  Subsection  6.2.1,  convective  vertical  velocity  is 
temporally  low-pass  but  has  a  non-zero  spectral  peak  on  the  wavenumber  axis  which 
indicates  considerable  spatial  variation.  At  a  low  vehicle  speed,  convection’s  spatial 
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Figure  6.16:  A  sample  time  series  and  its  PSD  generated  by  a  second-order  AR  model 
for  AUV  speed  0.1  m/s. 

variation  cannot  be  well  sensed  in  the  AUV-recorded  time  series.  The  measured 
process  is  basically  still  low-pass,  similar  to  the  internal  wave  measurement.  This 
similarity  is  quantitatively  described  by  the  mingled  spectra  of  the  two  processes  as 
shown  in  Figure  6.5  at  low  AUV  speeds. 

It  is  noted  that  in  internal  wave  fields  or  convection  fields,  the  ocean’s  horizontal 
advection  often  exists  [47],  [2],  The  advective  speed  can  be  larger  than  0.1  m/s  [47], 
For  a  mooring,  the  effect  of  advection  current  makes  its  measurement  equivalent  to 
that  made  on  a  low-speed  AUV.  Then  based  on  the  mooring  measurement,  it  would 
also  be  difficult  to  classify  vertical  velocities  of  convection  and  internal  waves. 

At  higher  AUV  speeds  0.25  m/s  and  1  m/s,  convection’s  spatial  peak  on  the  v- 
axis  is  apparently  projected  onto  the  /-axis  of  the  corresponding  mingled  spectrum, 
as  displayed  by  Figure  6.4.  Thus  convection’s  spatial  feature  is  brought  to  light  in 
the  AUV-recorded  time  series.  Due  to  the  properties  of  internal  wave’s  frequency- 
wavenumber  spectrum,  its  mingled  spectrum  lies  at  baseband  on  the  /-axis.  A  higher 
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Figure  6.17:  At  AUV  speed  1  m/s,  the  histograms  of  feature  z  for  internal  wave  and 
convection.  The  value  of  2  corresponding  to  the  experimental  data  is  marked  by  the 
arrow. 

AUV  speed  does  not  change  this  basic  spectral  shape.  Consequently,  the  mingled 
spectra  of  the  two  processes  show  a  noticeable  difference.  This  highlighted  difference 
greatly  improves  the  classifier’s  performance,  as  shown  by  the  much  better  separated 
histograms  and  the  corresponding  ROC  curves  in  Figure  6.18.  Fundamentally,  a 
higher  AUV  speed  pulls  the  peak  of  convection’s  mingled  spectrum  farther  away 
from  the  base  frequency  band  where  internal  wave’s  stays  despite  the  heightened 
vehicle  speed.  This  case  study  demonstrates  that  we  can  take  advantage  of  AUV’s 
controllable  motion  for  good  classification. 


Classification  feature  z 
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Table  6.2:  Statistics  of  z  (class  1:  internal  wave,  class  2:  convection) 


1  m/s 

mz  i 

°z\ 

mz2 

002 

Theoretical 

wm 1 

bniiiiE 

3.95 

Simulation 

2.42 

0.25  m/s 

mzi 

00l 

mz  2 

002 

Theoretical 

0.49 

1.24 

3.24 

1.99 

Simulation 

0.82 

0.86 

3.09 

1.84 

0.1  m/s 

mz  i 

00 1 

mz2 

0z2 

Theoretical 

ESI 

0.48 

0.42 

0.46 

Simulation 

0.27 

0.39 

0.36 

0.05  m/s 

mz  i 

00  1 

mz2 

002 

Theoretical 

-0.07 

iUEEl 

0.01 

0.21 

Simulation 

-0.06 

0.03 

0.23 
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Figure  6.18:  The  classifier’s  performance  at  a  series  of  AUV  speeds. 
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Chapter  7 


Classifier  Robustness 


7.1  Spectrum  Uncertainty  at  Successive  Stages 


Global  class  mean 

Myo 

Variant 

process  parameters  \ 

i 

r 

E[My]=Myo 

Local  class  mean 
My 

Uncertainty  from  ► 

unreliably  modeled  y 

r 

'  E[M\ 

',IMy]=My 

portion  ot  FSU 

Uncertain 
local  class  mean 
My’ 

Variance  in  . 

periodogram  estimate  1 

r 

E[SyIMy’]=My’ 

PSD  estimate 

Sy 

Figure  7.1:  Flow  chart  of  spectrum  uncertainties  at  successive  stages. 
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Essentially,  the  AUV-based  classifier  works  by  matching  the  input  data’s  spectrum 
with  the  model-based  spectrum  templates.  In  the  thesis,  the  two  models  are  the 
MIT  Convection  Model  and  the  Garrett-Munk  internal  wave  model,  as  introduced 
in  Chapter  5.  To  build  the  spectrum  templates,  we  assign  a  set  of  parameters  to  each 
model.  These  parameters  are  selected  based  on  our  understanding  of  the  process 
and  available  prior  information.  For  instance,  we  set  the  heat  flux  Q=300  W/m2 
and  mixed  layer  depth=350  m  for  the  convection  model,  based  on  meteorological  and 
hydrographic  measurements  in  the  Labrador  Sea. 

The  real  data,  however,  may  have  some  discrepancy  from  the  model.  From  the 
perspective  of  modeling,  this  mismatch  problem  is  equivalent  to  parameter  uncer¬ 
tainty.  We  should  enable  the  AUV-based  classifier  to  be  robust  to  model  parameter 
uncertainty.  As  briefly  mentioned  in  Section  6.3,  the  computation  of  the  feature 
transformation  vector  V  and  the  classification  results  presented  in  Section  6.5  have 
already  used  the  robust  classifier  presented  in  this  chapter.  We  did  not  present  ro¬ 
bust  design  therein  so  as  not  to  complicate  the  mainstream.  Adding  robustness  to 
the  AUV-based  classifier  becomes  the  topic  of  this  chapter. 

The  relation  between  the  model’s  template  spectrum  My0  and  the  data’s  PSD 
estimate  Sy  is  illustrated  in  Figure  7.1.  If  the  parameter  uncertainties  were  neglected, 
the  middle  two  blocks  would  not  exist  such  that  the  only  source  of  Var[Sy]  is  the 
periodogram’s  inherent  variance.  This  simplified  case  is  expressed  in  Equation  (4.23). 

Let  us  now  consider  parameter  uncertainties  associated  with  internal  waves.  As 
presented  in  Section  6.2,  the  buoyancy  frequency  N(z )  is  the  key  parameter  of  inter¬ 
nal  wave  vertical  velocity’s  spectrum  since  it  prescribes  the  upper  bound  frequency. 
While  we  set  this  parameter  to  about  3  cycles/hour,  we  should  allow  for  an  uncer¬ 
tainty  range  around  this  value.  Another  parameter  is  the  height  of  the  extended 

PSD  plateau  to  account  for  higher-frequency  processes.  It  is  set  such  that  the  ratio 

Power  above  buoyancy  frequency  i  n  i  i.  u  i  nr  ,  •  , 

- power  0f  internal  waves - -  equals  0.1.  We  should  also  allow  for  an  uncertainty 

range  for  this  ratio.  Using  internal  wave  as  an  example,  notations  in  Figure  7.1  are 

explained  as  follows. 

1.  My  and  Myo:  Under  a  specified  buoyancy  frequency  N(z),  the  mingled  spec- 
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trum  is  denoted  My  (including  the  data  window’s  effect).  Uncertain  values 
of  N(z),  however,  will  cause  major  changes  to  the  mingled  spectrum.  Corre¬ 
sponding  to  some  range  of  N(z ),  e.g.,  from  1  cycle/hour  to  5  cycles/hour,  MY 
experiences  a  variational  range  too.  The  mean  of  My  is  denoted  MY0: 


E[My]  —  My0 


(7.1) 


We  call  MY 0  the  “global  class  mean” . 

2.  My :  When  the  buoyancy  frequency  N(z)  is  specified,  so  will  be  My.  But 
spectrum  uncertainty  still  exists,  which  is  caused  by  the  contribution  from  the 
unreliably  modeled  portion  of  the  t]  —  v  PSD,  i.e.,  the  uncertain  height  of  the 
extended  PSD  plateau  at  above  the  buoyancy  frequency.  At  a  specified  plateau 
height,  the  mingled  spectrum  is  denoted  My.  Corresponding  to  some  range  of 
plateau  height  variation,  My  has  a  variational  range  too.  The  conditional  mean 
of  My  is  MY: 


E\My\My]  =  My 


(7.2) 


We  call  My  the  “local  class  mean” . 

3.  SY :  When  both  the  buoyancy  frequency  N(z )  and  the  height  of  the  extended 
plateau  are  specified,  My  will  be  specified.  Its  estimate  SY,  however,  has  un¬ 
certainty  due  to  the  nature  of  periodogram  estimation.  The  conditional  mean 
of  SY  is  Myl 


E[Sy\My]  =  My 


(7.3) 
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For  convective  vertical  velocity,  the  concept  is  the  same.  The  picture  is  simpler  by 
leaving  out  the  third  block  “uncertain  local  class  mean”:  as  we  regard  convection’s 
rj  —  v  PSD  as  reliably  modeled,  My  and  My  coincide,  or  equivalently,  V ar[My]  =  0. 

7.1.1  Local  Uncertainty 

As  presented  in  Section  6.2,  a  power  density  plateau  is  added  to  internal  wave  verti¬ 
cal  velocity’s  temporal-spatial  PSD  to  account  for  processes  at  above  the  buoyancy 
frequency.  By  calculations  based  on  a  notional  plot,  we  set  the  plateau  height  such 
that  its  integrated  power  equals  10%  of  internal  wave’s,  as  shown  in  the  lower  panel 
of  Figure  6.2.  In  contrast  to  internal  wave  PSD,  we  deem  the  extended  plateau  “unre¬ 
liably  modeled” .  It  contributes  to  the  mingled  spectrum  by  integration  as  illustrated 
in  Figure  3.2.  The  proportion  of  “unreliable  contribution”  varies  with  frequency  as 
the  integration  line  slides  through  the  temporal-spatial  PSD.  Figure  7.2  shows  the 
variation  of  percentage  of  internal  wave  power  density  contribution  (considered  reli¬ 
able)  as  a  function  of  frequency  and  AUV  speed.  At  a  lower  AUV  speed,  the  “reliable 
contribution  percentage”  drops  faster  to  zero  because  the  integration  line  in  Figure  3.2 
slides  out  of  the  reliable  PSD  region  earlier. 

There  exists  uncertainty  with  the  plateau  height,  which  implies  that  the  total 
power  at  above  the  buoyancy  frequency  could  vary.  This  uncertainty  translates  into 
the  mingled  spectrum  after  integration.  It  affects  the  PSD  estimate  variance  which 
in  turn  has  an  effect  on  the  transformation  vector  V . 

In  this  section,  we  do  not  yet  include  global  uncertainty,  but  only  consider  local 
uncertainty  due  to  the  variation  of  the  extended  plateau  height.  Following  notations 
in  Figure  7.1,  we  denote  internal  wave’s  uncertain  PSD  template  as  My,  and  its  mean 
as  My.  My  is  composed  of  reliable  and  unreliable  contributions: 


My  =  aMy  +  (1  -  a)My  (3  (7.4) 

where  a  is  the  reliable  contribution  percentage  as  shown  in  Figure  7.2.  Random 
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AUV  speed  =  0.25  m/s 


Figure  7.2:  Percentage  of  reliably  modeled  power  density  contribution  to  the  internal 
wave  mingled  spectrum. 

variable  /3  embodies  the  uncertainty  of  the  unreliable  contribution.  It  is  supposed  to 
have  a  mean  of  one: 


m  =  i 


(7.5) 


Then  we  have 


E[M y]  =  aMY  +  (1  -  a)MYE[p]  =  My 


(7.6) 


Var[My ]  =  My(  1  -  a)2Var[0\ 


(7.7) 
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These  local  uncertainty  properties  will  be  used  in  the  derivations  in  Section  7.2. 

7.1.2  Global  Uncertainty 

1.  Ocean  convection. 

In  Section  5.1,  we  set  the  heat  flux  Q=300  W/m2  and  the  mixed  layer  depth 
H= 350  m  for  the  convection  model,  based  on  meteorological  and  hydrographic 
measurements  in  the  Labrador  Sea.  Corresponding  to  the  concept  in  Figure  7.1, 
we  regard  the  mingled  spectrum  under  those  parameters  as  the  global  class 
mean  My2_o-  Now  we  apply  a  significantly  different  set  of  parameters  to  repre¬ 
sent  global  uncertainty.  We  let  the  heat  flux  and  the  mixed  layer  depth  both 
increase  by  a  factor  of  three:  Qextreme= 900  W/m2  and  //extreme— 1050  m.  We 
call  convection  under  this  setting  an  “extreme  convection  case” ,  and  denote  the 
resultant  mingled  spectrum  Myi-extreme- 

At  36000  s  (10  hours)  after  surface  cooling  starts,  vertical  flow  velocity  w  is 
shown  in  Figure  7.3.  The  upper  panel  displays  the  horizontal  cross-section  at  the 
500-m  depth;  the  lower  panel  displays  the  vertical  cross-section  at  y  —  1000  m. 
Convective  cells  with  periodicity  of  about  400  m  ~  450  m  are  observable  in  both 
panels.  This  spatial  period  is  about  twice  the  value  (200  m  ~  250  m  as  shown 
in  Figure  5.2)  for  the  former  model  setting.  According  to  Equation  (5.2),  the 
increment  ratio  of  spatial  period  is  supposed  to  be  \/3  (noting  that  the  buoyancy 
flux  Bo  is  proportional  to  the  heat  flux).  The  new  model  run  result  is  close  to 
the  theoretical  prediction. 

Based  on  the  model  output  for  two  hours  (from  8  hours  to  10  hours  after  the 
onset  of  surface  cooling),  we  compute  the  temporal-spatial  PSD  of  convective 
vertical  velocity.  It  is  shown  in  the  upper  panel  of  Figure  7.4.  On  the  zA-axis, 
the  spectral  peak  wavenumber  is  lower  than  that  of  the  former  convection  result 
(shown  in  the  upper  panel  of  Figure  6.2),  due  to  an  enlarged  spatial  period  of 
convective  cells. 

2.  Internal  waves. 
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W  (m/s).  Time  step  3600.  Initially  uniform  profiles. Depth  500m. 


-0.5  -0.4  -0.3  -0.2  -0.1  0  0.1  0.2  0.3 


Y=1000m. 


Figure  7.3:  Horizontal  and  vertical  cross-sections  of  vertical  flow  velocity  w.  at.  time 
10  hours,  for  “extreme”  convection. 

In  Section  5.2,  we  set  N(z )  =  5.0  x  10-3  rad/s,  i.e.,  about  3  cycles/hour.  We 
regard  the  mingled  spectrum  under  this  parameter  as  the  global  class  mean 
Myi_o-  To  represent  global  parameter  uncertainty,  we  now  set  a  significantly 
different  value:  N(z )  =  1.7  x  10~3  rad/s,  i.e.,  about  1  cycle/hour.  We  call  in¬ 
ternal  wave  under  this  setting  an  “extreme  internal  wave  case” ,  and  denote  the 
resultant  mingled  spectrum  My1_ea;treme.  The  corresponding  temporal-spatial 
PSD  of  internal  wave  vertical  velocity  is  shown  in  the  lower  panel  of  Figure  7.5. 
Due  to  the  decreased  buoyancy  frequency,  the  internal  wave’s  spectrum  shrinks 
on  the  77-axis.  As  before,  we  have  made  an  extension  of  PSD  plateau  to  ac- 
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lOlog  S(ri,v)  for  convective  vertical  velocity 
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Figure  7.4:  Temporal-spatial  PSD  of  vertical  velocity  of  “extreme”  convection  (upper) 
and  “mean”  internal  wave  (lower). 

count  for  higher-frequency  processes.  The  plateau  height  is  set  such  that  the 

ratio  Power  above  buoyancy  frequency  ls  0  j  Note  that  the  internal  wave 
Power  of  internal  waves  '  -  . 

vertical  velocity’s  power  has  been  normalized  to  that  of  convection. 


-  From  the  rj  -  u  spectra,  we  compute  the  “mean”  and  “extreme”  mingled  spectra 
MYi_ o  and  MYi_extreme  (both  including  the  data,  window’s  effect).  Parameters  for 
MYi_a  and  MYi_extreme  are  summarized  in  Table  7.1.  At  AUV  speed  1  m/s,  the  two 
mingled  spectra  for  convective  vertical  velocity  are  shown  in  Figure  7.6,  along  with 
their  difference.  For  the  extreme  convection,  the  spacing  between  convection  cells  is 
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Figure  7.5:  Temporal-spatial  PSD  of  vertical  velocity  of  “mean”  convection  (upper) 
and  “extreme”  internal  wave  (lower). 

larger.  Consequently,  by  mingled  spectrum  projection  as  illustrated  in  Figure  6.4,  its 
spectral  peak  lies  at  a  lower  frequency. 

At  a  series  of  AUV  speeds  u=  1  m/s,  0.25  m/s,  0.1  m/s,  and  0.05  m/s,  the  “mean” 
and  “extreme”  mingled  spectra  as  well  as  their  difference  are  shown  in  Figure  7.7. 
Corresponding  to  the  second  block  in  Figure  7.1,  we  consider  (MYi_o  —  My]  extreme)2 
and  (MY2_ o  —  extreme)2  to  be  the  estimates  of  class  variances,  i.e.,  at  frequency 
k, 
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Table  7.1:  Parameters  for  MYl_o  and  My  Extreme 


Mean:  MYi_0 

Extreme.  ■^■Yi^e.xtreme 

Internal  wave  vertical  velocity 
(*=1) 

N(z )  =  5.0  x  10~3  rad/s 

N0  =  5.2  x  lO"3  rad/s 

N(z)  =  1.7  x  lO”3  rad/s 

N0  =  5.2  x  10“3  rad/s 

Convective  vertical  velocity 

(i= 2) 

surface  heat  flux=300  W/m2 
mixed  layer  depth=350  m 

surface  heat  flux=900  W/ m2 
mixed  layer  depth=1050  m 

Var[Myi(k)\  ft*  (MYuo(k)  —  Myi  -i extreme  (k)f  (7.8) 


V  ar[MY2{k)\  «  (MY2_0(k)  —  MY2  -extreme  (k))2  (7.9) 

7.2  Derivation  of  Total  Variance  of  Mingled  Spec¬ 
trum  Estimate 

With  considerations  of  local  and  global  model  parameter  uncertainties,  we  need  to 
re-examine  the  statistics  of  PSD  estimate  SY  of  the  AUV-recorded  time  series  Y(t). 
Only  on  this  basis  can  we  correspondingly  modify  the  feature  transformation  vector 
V  and  thus  add  robustness  to  the  classifier.  Note  that  the  following  derivations  are 
for  an  arbitrary  frequency  point. 

First,  let  us  find  the  expectation  of  SY: 
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At  AUV  speed  1  m/s 


Figure  7.6:  The  “mean”  and  “extreme”  mingled  spectra,  and  their  difference,  of 
convective  vertical  velocity  at  AUV  speed  1  m/s. 
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where  we  utilize  5y  P5K|M^(5v|My)  d5y  =  £’[5y|My]  =  My  (Equation  (7.3)), 


and  My  ,M  (My|My)  dMy  =  E[My|My]  —  My  (Equation  (7.2)),  as  well 
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Internal  wave.  AUV  speed  =  1  m/s  Convection.  AUV  speed  =  1  m/s 


AUV  speed  =  0.25  m/s  AUV  speed  =  0.25  m/s 


AUV  speed  =  0.05  m/s 


Frequency  (Hz) 


AUV  speed  =  0.05  m/s 


Frequency  (Hz) 


Figure  7.7:  The  “mean”  and  “extreme”  mingled  spectra  and  their  difference,  for 
internal  wave  and  convection  at  a  series  of  AUV  speeds. 


as  Equation  (7.1)  such  that  My  Pmy(^y)  dMy  —  Myo 
'  So  we  see  that  based  on  the  formulation  in  Section  7.1,  Sy  is  still  an  unbiased 
estimate  of  the  global  class  mean  Myo.  Next  let  us  find  the  variance  of  Sy.  Us¬ 
ing  Equation  (7.10),  we  have  , 


Var[SY]  =  E[(Sy  -  £[5y])2] 

=  £[(5y  -  Myo)2] 

(Sy  —  Myo)2  P§Y(Sy)  dSy 

[«y  —  M‘y )  +  (My  —  My)  +  {My  —  Myo)]2  Jlyy  (Sy)  dSy 

v2  (7-ID 
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[(5k  -  My)2  +  (My  -  My)2  +  (My  -  Myof 

+  2(My  -  My)(5y  -  M'y)  +  2 (My  -  My0)(5y  -  M^) 
+  2(My  -  My )  (-My  -  -My0)]  (5y )  dSy 


Let  us  integrate  the  six  terms  individually  and  then  sum  them  up.  Beginning  from 
the  integration  of  the  first  term: 
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where1 


r00  „  ,  vr'2 

/  (5y  -  M(y)2^  ,M,  (5y|M(y)  d5y  =  Var[Sy|M(,]  =  (7.13) 

J —oo  V eff 

Incorporating  Equation  (7.13)  into  Equation  (7.12),  we  have 


1The  zero-frequency  is  special:  V  ar[SY  (S)\M'Y]  —  2--f^ .  This  specialty  is  properly  treated  in 
computations,  but  omitted  in  expressions  herein  for  the  sake  of  conciseness. 


(7.14) 


—  My)2  Psy(Sy)  dSy 
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OO 

/°°  (My  -  -E[My  |My]  +  £[M(,|My])2  PM^|Ary  (M^|My) 

J  — OO 

J  |  (2?[My|My])2  +  (My  —  S[My|My])2 
+  2J5[My|My](My  -  E[M'y\My})  |  PM^|MK(M^|My)  dMy 


(7.15) 


By  Equation  (7.2),  the  first  term  equals  My;  the  second  term  equals  V  ar[My\My]\ 
the  third  term  vanishes.  Note  that  in  Equation  (7.2),  “. . .  |My”  was  implicitly  in¬ 
cluded  as  My  was  still  considered  deterministic  therein.  Equation  (7.15)  then  becomes 


My  PM^jMy  (My  |My)  dMy 


=  My  +  Var[My\My] 

=  My  +  (l  -  a(MY))2Var[f3\M$- 
=  {1  +  (1  -  a(My))Var[/3]}Mj 


(7.16) 


where  we  utilize  Equation  (7.7)  to  reduce  Var[My\My].  Once  again,  note  that 
in  Equation  (7.7),  “. . .  |My”  was  implicitly  included.  In  Equation  (7.16)  we  assume 
/?  does  not  vary  with  My  but  a  does,  thus  we  use  notation  a(My).  Then  Equa¬ 
tion  (7.14)  is  written  as 


123 


1 


0° 

{l  +  (l  —  Oi(My))  Var[(3]}My  Pmy(My)  dMy 

UeffJ- oo 


a  (percentage  of  reliably  modeled  PSD  contribution  for  Internal  wave).  AUV  speed  =  1  m/s 


Frequency  (Hz) 


Figure  7.8:  Internal  wave’s  a  for  the  “mean”  and  “extreme”  cases. 


For  the  “mean”  and  “extreme”  internal  wave  cases,  a  (percentage  of  the  reliably 
modeled  PSD  contribution)  is  shown  in  Figure  7.8.  a  drops  with  frequency  as  the 
rningled  spectrum  integration  line  slides  out  of  the  central  internal  wave  spectrum, 
into  the  outside  plateau.  The  rate  of  descent,  however,  varies  with  the  AUV  speed 
since  the  integration  line’s  slope  is  determined  by  the  vehicle  speed.  At  a  given  AUV 
speed,  a  is  lower  for  the  “extreme”  case  than  for  the  “mean”  case  because  the  central 
internal  wave  spectrum  is  narrower  in  the  former.  Variational  a(MY)  adds  complexity 
to  the  integration  in  Equation  (7.17).  For  the  sake  of  robustness  for  classification, 
let  us  find  an  upper  bound  of  the  integral.  We  denote  the  minimum  value  of  a 


corresponding  to  the  range  of  My  as  amin.  At  any  vehicle  speed,  we  pick  the  lower 
one  of  the  two  a  curves  (in  Figure  7.8)  as  amin.  Then  we  have 


1  +  (l  -  a(MY))2Var[P ]  ^  1  +  (1  -  amin)2Var[j3 ]  (7.18) 
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/: 


My  pmy  [My  )  dMY 


where 


My  pmy  {My)  dMy 


=  ( E[MY ])2  +  Var[MY } 
=  M2  o  +  Var[MY ] 


(7.20) 


where  we  utilize  Equation  (7.1)  such  that  E[MY\  =  MY0. 

Combining  Equation  (7.17),  Equation  (7.19),  and  Equation  (7.20),  we  reduce  Equa¬ 
tion  (7.14)  to 


r  (Sr  -M'yf  Ph  (Sy)  dSY  S  {l  +  (1  a^Var[^  (Mj0  +  Var[Mr\) 

J  — oo  ^eff 

(7.21) 

We  have  thus  obtained  an  upper  bound  of  the  integration  of  the  first  term  in  Equa¬ 
tion  (7.11).  Now  let  us  find  the  integration  of  the  second  term: 
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where  /^0Pj5y|M^(5y|My)dS'y  =  1.  So  Equation  (7.22)  is  simplified  to 
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=  (1  -  <w)Var[/?]((£[My])2  +  Var[MY ]) 

=  (1  -  amin)2Var[P](M20  +  Ear[My]) 


(7.23) 


where  we  utilize  E[Miy\MY ]  =  My  (Equation  (7.2))  and  Ear[M(,|My]  =  (l  - 
a(My))  Ear[/?]My  (Equation  (7.7)).  The  same  as  in  the  derivation  of  Equation  (7.21), 
we  use  a  ^  amin  to  arrive  at  an  upper  bound  of  the  integral. 

Similarly,  we  can  find  the  integration  of  the  third  term  in  Equation  (7.11): 
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where  we  note  that  JZoPsy\m^(Sy\My)  d$Y  =  1,  !T00Pm^\my(M'y\My)  dM’y  =  1, 
and  that  E[MY]  —  MY0. 

Now  let  us  show  that  the  integrations  of  the  fourth,  fifth,  and  sixth  terms  in  Equa¬ 
tion  (7.11)  all  vanish.  Integration  of  the  fourth  term  is 
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since  we  know  that  f^X)( Sy  —  My)  p m<y j My  (M'Y\ My )  dSy 
=  E[(Sy  -  My)\My\  =M!y~M'y=  0. 

Integration  of  the  fifth  term  in  Equation  (7.11)  is 
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/  2{My  —  Myo)(SY  —  My)  PsY(Sy)  dSy 

J  —  oo 

/oo 

{My  -  My0)(Sy  -  M'y) 

-oo 

roc  ^  poo 

{  PsYW^SylM'y)  dM'y  /  PM'y\My  (My\My)  pMy  (My)  dMy}  dSy 
j  -oo  J  —  00 

/oo  poo 

{My  —  Myo)  Pmy(MY)  dMy  /  PM'Y\MY{My\My)  dM'y  2gj 

/oo 

(S'y  —  My)  P§y\M'(Sy\My)  dSy 

•oo 

=  0 


due  to  the  the  same  reason  as  for  Equation  (7.25). 
Integration  of  the  sixth  term  in  Equation  (7.11)  is 


/OO 

2 (My  —  My )  (My  —  Myo)  P<jy  (Sy)  dSy 

-oo 

/oo 

(My  -  My) (My  -  Myo) 

-oo 

r°°  ^  r°° 

{  /  P.Sy|M^  (‘S’y  |My)  dMy  /  I  My  (My  |My )  pMy  (My)  dMy}  dSy 

J  —  oo  J  — OO 

/oo 

(My  —  Myo)  PMy  (My)  dMy 

-oo 

r°°  r°o  (7.27) 

/  (M!y-My)pM,Y\MY(M'Y\MY)dM'y  I  PsyW>  (Sy\M!y)  dSy 

•/-oo  d-oo  y 

/°°  /‘OO 

(My  -  Myo)  PMy  (My)  dMy  /  (My  -  My)  Pm'  |My  (My  |My)  dMy 
oo  J  —  oo 


=  0 


noting  that  p^|M^  (5y  |M^)d5y  =  1  and  that  (M|, -My  )pM,.,My  (M^My)dM', 

E’[(My  —  My)|My]  =  My  —  My  =  0 

Summing  up  Equation  (7.21),  Equation  (7.23),  Equation  (7.24),  Equation  (7.25), 
Equation  (7.26),  and  Equation  (7.27),  Equation  (7.11)  finally  reduces  to 
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Var[SY]  <  {(1  -  aminfVar[P]  +  1  +  t1  <^?Var\fi\  ^ 

Veff 

+  (1  +  ){l  +  (1  -  amin)2V  ar[/3]}V  ar[My] 

Veff 


(7.28) 


Equation  (7.28)  expresses  the  total  variance  of  the  PSD  estimate  Sy ,  including 
effects  of  local  uncertainty  (embodied  in  (1  —  amin)2Par[^])  and  global  uncertainty 
(embodied  in  V ar[My])  of  model  parameters.  Let  us  observe  Equation  (7.28)  in 
simplified  cases: 


1.  When  not  considering  global  variation  of  My,  i.e.,  Var[My]  =  0  and  thus  a  has 
no  variation,  Equation  (7.28)  is  simplified  to 


Var[SY]  =  {(1  -  a)2Var[p]  +  1  +  ^ . ^~V(lr^}M^0 

Ueff 

(7.29) 

This  corresponds  to  cases  where  only  local  uncertainty  is  a  concern. 

2.  When  a  =  1  or  Var[/3]  =  0,  i.e.,  either  the  unreliably  modeled  PSD  makes  no 
contribution  or  the  PSD  is  reliably  modeled,  local  uncertainty  is  not  a  concern 
(as  for  convection).  Then  Equation  (7.28)  is  simplified  to 


Var[SY }  =  —  M20  +  (1  +  -1 ~)Var[MY }  (7.30) 

ueff  Veff 

3.  If  not  only  local  uncertainty  is  left  out  but  also  we  have  veff  1  ,  then  Equa¬ 
tion  (7.30)  is  further  simplified  to 


Var[Sy\  Var[My] 


(7.31) 
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i.e.,  only  V ar[MY]  leaves  its  legacy. 


4.  When  both  local  and  global  uncertainties  are  neglected,  Equation  (7.30)  is  sim¬ 
plified  to 

Var[SY ]  =  —  M£0  (7.32) 

Ueff 

This  is  the  most  simplified  form  of  Equation  (7.28),  only  containing  the  pe- 
riodogram’s  inherent  variance.  We  have  seen  it  in  Equation  (4.23).  Compar¬ 
ing  Equation  (7.32)  with  Equation  (7.28),  we  note  that  the  model  parameter 
uncertainties  increase  Var[SY ]  in  addition  to  the  periodogram’s  inherent  vari¬ 
ance. 

7.3  Modified  Feature  Transformation  Vector 

With  an  increased  Var[SYi(k)],  the  formulation  of  the  transformation  vector  V  in  Sub¬ 
section  4.3.1  needs  to  be  modified  in  the  steps  as  follows. 

1.  Modify  covariance  matrices  Eyi  and  Ey2. 

The  formulation  of  the  covariance  matrix  of  SYi  in  Equation  (4.22)  still  holds, 
but  we  should  assign  updated  values  to  the  diagonal  elements.  For  internal 
wave  vertical  velocity,  the  covariance  matrix  is 


Eyi  =  diag{Var[SY1( 0)],  •  •  •  ,  Ear[Syi(iV  -  1)]} 

(7.33) 


where  (0),  •  •  •  ,  (TV  — 1)  denote  the  N  frequency  points.  At  frequency  k,  diagonal 
element  Var[5yi(fc)]  is  replaced  by  its  upper  bound  as  in  Equation  (7.28): 
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Var[Sri(k)}  =  {(1  -  ami„(k))2Var{0}  +  1  +  t1  „(*) 

Ueff 

+  (1  +  ~~){1  +  C1  “  otmin{k))2V ar[(3\}V ar[MYi]  ^  ^ 

Similarly,  for  vertical  velocity  of  convection,  the  covariance  matrix  is 


Ey2  =  diag{Var[SY2( 0)],  •  •  -  ,  Var[SY2(N  -  1)]} 

(7.35) 

At  frequency  k,  diagonal  element  Var[SY2(k)}  is  replaced  by  its  upper  bound 
as  in  Equation  (7.28): 


Var[SY2(k)}  =  Q(k)  +  (1  +  — )Var[MY2 1  (7.36) 

Veff  Veff 

noting  that  we  regard  convection’s  77  —  u  PSD  as  reliably  modeled  (i.e.,  there 
exists  no  local  uncertainty)  such  that  amin  =  1,  which  simplifies  the  expression. 

Mingled  spectrum  class  means  MY i_0  (for  internal  wave  vertical  velocity)  and 
My 2_o  (for  convective  vertical  velocity)  are  shown  by  the  blue  curves  in  Fig¬ 
ure  7.7.  For  the  extreme  internal  wave  and  convection  cases,  MYi_extreme  and 
MY2 -extreme  are  shown  by  the  red  curves  in  Figure  7.7.  V ar[MY  1]  and  V ar\MY 2] 
are  then  calculated  by  Equation  (7.8)  and  Equation  (7.9),  based  on  the  spec¬ 
trum  difference  shown  by  the  crosses  in  Figure  7.7. 

Knowing  convection’s  MY2_ 0  and  Var[MY2]  is  sufficient  for  computing  Var[SY2] 
by  Equation  (7.36).  To  compute  internal  wave’s  Var[SY\ ]  by  Equation  (7.34), 
however,  we  should  also  know  otmin  and  Var[(3\  (as  in  Equation  (7.4))  besides 
My i_0  and  Var[MYi].  a  (percentage  of  the  reliably  modeled  PSD  contribution) 
for  the  “mean”  and  “extreme”  internal  wave  cases  is  shown  in  Figure  7.8.  At  any 
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vehicle  speed,  we  pick  the  lower  one  of  the  two  a  curves  as  amin.  The  extended 
plateau  PSD  height’s  variance,  Var[(3],  is  set  to  4.  This  value  is  considered 
sufficiently  large  to  account  for  the  plateau  height’s  uncertainty  (please  refer 
to  Subsection  6.2.1  for  discussions).  Now,  with  My i_0,  Var[MY1],  amin,  and 
Var[(3 ]  all  known,  we  can  compute  internal  wave’s  Uar[Syi]  by  Equation  (7.34). 

2.  Update  transformation  vector  V. 

With  Uar[5yi]  and  V  ar[SY2]  known,  Eyi  and  Ey2  are  obtained  by  Equa¬ 
tion  (7.33)  and  Equation  (7.35),  respectively.  Incorporating  E y\  and  Ey2 
into  Equation  (4.4),  we  get  the  within-class  scatter  matrix  Aw_y  (assuming 
Pi  =  P2  =  f).  Then,  incorporating  Aw_y,  MY1_0,  and  My i_0  into  Equa¬ 
tion  (4.2),  we  arrive  at  the  transformation  vector  V: 


V  —  Aw\y  (My 2_o  —  Myx_ 0) 

_  My  2_o(0)  —  Myi_o(0) 

~  -5[yar[^i(0)]  +  ^fe(0)]]  (7-37) 

My2_o(iV  —  1)  —  MY1_o(N  —  1)  T 
\  \y o-'f[SYi(N  -  1)]  +  Var[SY2(N  -  1)]] . 

After  including  model  parameter  uncertainties,  the  total  variance  of  SYi  is  in¬ 
creased  compared  with  merely  the  periodogram’s  inherent  variance.  Conse¬ 
quently,  the  denominators  in  Equation  (7.37)  are  larger  than  their  counterparts 
in  Equation  (4.25).  The  magnitude  of  the  modified  V  is  therefore  smaller  than 
that  of  the  unmodified  V,  as  illustrated  in  Figure  7.9.  For  AUV  speed  1  m/s, 
the  modified  and  the  unmodified  transformation  vectors  are  compared  in  the 
lower  panel  of  Figure  6.6.  A  lower  magnitude  of  V  means  that  to  a  lesser  extent 
we  utilize  the  difference  between  the  two  spectra.  The  classifier’s  performance 
is  thus  lowered,  but  with  a  gain  of  robustness. 
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Figure  7.9:  The  impact  of  a  larger  spectrum  variance  on  the  transformation  vector 
formulation. 

7.4  Robustness  Tests  on  AUV-Based  Classifier 


We  test  the  classifier’s  robustness  using  input  data  that  is  mismatched  with  the 
model.  The  goal  is  to  distinguish  the  “extreme”  convection  from  the  “mean”  internal 
wave.  As  the  “mean”  internal  wave  spectrum  has  a  higher  buoyancy  frequency  (about 
3  cycles/hour)  than  that  of  the  “extreme”  one  (about  1  cycle/hour),  its  mingled  spec¬ 
trum  is  closer  to  that  of  convection.  Hence  the  classification  is  more  challenging.  We 
therefore  select  this  more  difficult  case  to  test  the  classifier’s  robustness.  Test  case 
parameters  and  approaches  are  summarized  in  Table  7.2,  as  compared  against  Ta¬ 
ble  6.1. 


Table  7.2:  Parameters  and  method  for  classifier  test  of  the  “extreme”  convection  case 
versus  the  “mean”  internal  wave  case 


Test  parameters 

Simulation  method 

Internal  wave  vertical  velocity 
(mean  case) 

N(z)  =  5.0  x  10~3  rad/s 

Nq  =  5.2  x  10~3  rad/s 

200  time  series  by 

AR  modeling 

Convective  vertical  velocity 
(extreme  case) 

200  time  series  by 

extracting  AUV  survey  lines  from 

the  convective  box  at  depth  500  m 
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1.  I Jeff  ~  1. 

When  the  AUV-acquired  data  is  short,  it  will  be  hard  to  do  time-domain  seg¬ 
mentation  or  frequency-domain  smoothing  for  the  purpose  of  classification.  This 
situation  necessitates  letting  ueff  =  1  in  formulating  the  classifier.  Detailed  dis¬ 
cussions  are  found  at  the  end  of  Section  6.3. 

The  modified  transformation  vectors  for  a  series  of  AUV  speeds  have  been  shown 
in  Figure  6.7,  and  explained  in  Section  6.3.  The  classifier’s  performance  is 
shown  in  Figure  7.10.  We  see  that  with  the  modification  of  V,  the  AUV-based 
classifier  is  able  to  distinguish  model-mismatched  convection  from  internal  waves 
at  a  high  AUV  speed.  As  we  have  observed  before,  the  classifier’s  performance 
improves  with  the  AUV  speed,  since  a  higher  vehicle  speed  pulls  the  peak  of 
convection’s  mingled  spectrum  farther  away  from  the  base  frequency  band  where 
internal  wave’s  mingled  spectrum  stays. 

2.  veff  —  4. 

If  a  longer  data  record  is  available,  we  can  apply  time-domain  segmentation  or 
frequency-domain  smoothing  to  increase  ueff  so  as  to  reduce  the  periodogram’s 
variance.  Now  we  assume  that  the  data  record  is  lengthened  to  5600  s=1400  s 
x  4.  We  then  partition  the  time  series  into  four  non-overlapping  segments,  and 
use  the  average  of  the  four  periodograms  as  the  input  to  the  classifier.  Thus 
Veff  equals  4,  which  reduces  the  periodogram’s  variance  by  a  factor  of  4. 

A  change  of  vejj  updates  Var [5Vi]  and  Vor[5,y2]  according  to  Equation  (7.34) 
and  Equation  (7.36),  respectively.  At  AUV  speed  1  m/s,  the  “modified”  V  is 
shown  in  the  lower  panel  of  Figure  7.11.  Compared  with  the  “unmodified”  V 
(i.e.,  only  taking  into  account  the  periodogram’s  inherent  variance  but  without 
consideration  of  the  spectrum’s  local  or  global  uncertainty),  the  “modified”  V 
has  a  lower  magnitude  because  of  the  additional  variances.  Compared  with  Fig¬ 
ure  6.6  where  ueff  =  1,  the  magnitudes  of  both  the  “unmodified”  V  and  the 
“modified  V”  are  higher  due  to  the  larger  vefj.  With  reduction  of  the  peri¬ 
odogram’s  inherent  variance,  variances  due  to  the  spectrum’s  local  uncertainty 
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(only  for  internal  wave  because  of  the  unreliably  modeled  PSD  plateau)  and 
global  uncertainty  (for  both  internal  wave  and  convection  because  of  signifi¬ 
cant  variations  of  process  parameters)  carry  more  weight  in  the  total  variance. 
Consequently,  the  effect  that  variances  from  the  spectrum’s  local  and  global  un¬ 
certainties  restrain  the  growth  of  V  is  more  apparent  in  Figure  7.11  (veff  =  4) 
than  in  Figure  6.6  (veff  =  1)- 

Corresponding  to  Figure  6.7,  V  for  a  series  of  AUV  speeds  is  shown  in  Fig¬ 
ure  7.12  when  veff  =  4.  Then  we  test  the  classifier  again  by  the  “extreme” 
convection  and  the  “mean”  internal  wave  cases.  Test  case  parameters  and  ap¬ 
proaches  are  the  same2  as  in  Table  7.2.  The  difference  is  that  we  now  take 
time  series  four  times  as  long,  and  then  use  the  averaged  periodogram  of  the 
four  segments  as  the  input  to  the  classifier.  The  classifier’s  performance  is 
shown  in  Figure  7.13.  We  see  that  compared  with  vejf  =  1  in  Figure  7.10, 
the  classifier’s  performance  improves  at  AUV  speed  1  m/s  and  0.25  m/s.  This 
improvement  is  due  to  the  reduced  variance  of  the  periodogram  because  of  the 
increased  ueff.  At  low  AUV  speeds  of  0.1  m/s  and  0.05  m/s,  however,  the 
classifier  still  has  a  hard  time  distinguishing  the  two  processes.  This  is  because 
the  two  mingled  spectra  are  similar  at  low  AUV  speeds.  Simply  reducing  the 
periodogram’s  variance  cannot  essentially  break  this  similarity.  We  still  need  to 
resort  to  a  higher  vehicle  speed  for  a  better  classification. 


2Except  that  at  AUV  speed  1  m/s,  100  (rather  than  200)  time  series  are  tested  because  of  the 
limitation  imposed  by  the  convective  box’s  size. 
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1  m/s 


0.25  m/s 


0.1  m/s 


Class  2:  convactlpn.  AUV  speed=0.25  m/s.  m  =2.91  .  a  =1.0 


Classification  feature  z 

Class  1:  internal  wave.  AUV  speed=0. 1  m/s.  m  =0.23,  o  =0.31 


Classification  feature  z 

Class  2:  convection.  AUV  speed=O.OS  m/s.  m^=-0.04,  <j  =0.28 


Classification  feature  z 


Figure  7.10:  The  classifier’s  performance  for  the  “extreme”  convection  versus  the 
“mean”  internal  wave,  when  ueff  —  1. 
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AUV  speed  =  1  m/s 


Frequency  {Hz) 


Frequency  (Hz) 


Figure  7.11:  At  AUV  speed  1  m/s,  class  mean  vectors  MY\_ o  and  My2_o  (upper 
panel),  and  the  transformation  vector  V  (lower  panel)  when  ueff  =  4. 


'  (when  ueff  =  4)  at  a  series  of  AUV  speeds. 
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Chapter  8 


AUV-Borne  Flow  Velocity 
Measurement  and  Data  Processing 

8.1  System  Integration 

An  Acoustic  Doppler  Velocimeter  (ADV)  [68]  measures  water  current  velocity  utilizing 
the  Doppler  principle  [35].  It  transmits  acoustic  waves  and  then  receives  echoes 
returning  from  sound  scatterers  in  the  water.  The  reflected  wave  has  a  frequency 
shift  compared  with  the  transmitted  wave.  The  frequency  shift  is  proportional  to  the 
radial  velocity  of  the  scatterer,  as  expressed  by  Equation  (8.1). 


Id  —  2 —fs  (8.1) 

c 

where  fr>  is  the  Doppler  frequency  shift;  fs  is  the  frequency  of  the  transmitted  signal; 
vr  is  the  radial  velocity  of  the  scatterer;  c  is  the  sound  speed.  Note  that  since  an 
ADV  both  transmits  and  receives,  the  Doppler  frequency  shift  is  doubled. 

Based  on  the  frequency  difference  between  the  transmitted  and  received  signals, 
the  velocity  of  sound  scatterers  can  be  calculated.  These  scatterers  are  plankton 
or  other  small  particles  floating  in  the  water.  In  most  cases,  the  assumption  that 
the  scatterers  are  passively  advected  by  water  motion  is  valid,  hence  the  scatterers’ 
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velocity  represents  that  of  the  water  current  [69],  [70]. 


Acoustic  receivers 


Acoustic  transmitter 


Figure  8.1:  Side  view  of  an  Acoustic  Doppler  Velocimeter  (ADV)  probe  (based  on 
Figure  4  of  [4]). 

An  ADV  probe  is  illustrated  in  Figure  8.1.  The  acoustic  beams  of  the  transmitter 
and  the  three  receivers  intersect  at  a  small  sampling  volume  (<  2  cm3)  located  away 
from  the  instrument  base  (16  cm  distance  for  Model  ADVOcean.  Please  see  Table  8.1). 
Three-dimensional  flow  velocity  at  this  distant  focal  point  is  calculated.  Thus  the 
flow  velocity  measurement  can  be  considered  undisturbed  by  the  probe.  As  shown 
in  Figure  8.1,  the  ADV’s  local  z-axis  is  defined  along  the  probe  stem;  the  z-axis  is 
coplanar  with  one  designated  receiver  arm;  the  y-axis  is  accordingly  defined  by  the 
right-hand  rule.  Table  8.1  gives  the  specifications  [4]  of  the  ADV  device  that  we  have 
installed  in  an  Odyssey  IIB  AUV. 

An  ADV’s  spatial  focus  and  low  noise  make  it  suitable  for  experiments  that  require 
high-resolution  and  high-precision  [71].  It  has  found  applications  in  flow  measurement 
in  laboratory  flumes  [72],  near  river  beds  [73],  [74],  and  the  seabed  [75].  In  the  above 
applications,  ADVs  monitor  current  velocity  only  at  spatially  fixed  positions.  The 
1998  Labrador  Sea  Experiment  requires  high-precision  flow  measurement  of  weak 
velocity  signals  (several  cm/s)  from  an  AUV.  This  provides  a  unique  opportunity  and 
challenge  of  integrating  an  ADV  into  a  moving  platform. 

The  mounting  location  and  orientation  of  an  ADV  on  an  AUV  should  meet  the 
following  requirements: 

1.  Avoid  saturation  of  velocity  measurement. 
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Table  8.1:  Specifications  for  the  SonTek  ADVOcean 


Sound  frequency 

5  MHz 

Output  data  rate 

0.1  Hz  ~  25  Hz 

Velocity’s  dynamic  range  (x  and  y) 

0  to  0.5,  1.2,  2.0,  6.0,  7.2  m/s,  programmable. 

Velocity’s  dynamic  range  (z) 

1/4  of  above 

Measurement  noise 
(at  25  Hz  data  rate) 

1%  of  velocity  range 

Velocity  resolution 

10-4  m/s 

Distance  of  sampling  volume  from  transmitter 

0.16  m 

Sampling  volume  size 

2  cm3 

Depth  rating 

2000  m 

Size 

0.36  m  x  0.18  m  (diameter  of  stem:  0.05  m) 

Weight 

1.5  kg 

2.  Avoid  interference  with  other  AUV  instruments. 

3.  Keep  the  ADV  probe  out  of  harm’s  way  during  AUV  launch  and  recovery. 

4.  Minimize  hull  influence  corrections  so  that  the  measurements  are  as  direct  as 
possible. 

5.  Locate  the  sampling  volume  outside  wakes  as  much  of  the  time  as  possible. 

According  to  the  above  requirements,  the  following  options  are  ruled  out:  1.  At  AUV’s 
nose.  The  vehicle’s  up  to  2  m/s  speed  could  saturate  the  ADV’s  z-velocity  with  the 
probe  in  the  along-ship  direction.  The  ADV  would  also  interfere  with  and  be  affected 
by  the  AUV’s  docking  latch  and  the  Ultra-Short-BaseLine  (USBL)  hydrophone  array 
that  are  mounted  at  the  nose.  The  velocity  measurement  would  need  considerable 
corrections  because  of  the  ADV’s  closeness  to  the  vehicle’s  stagnation  point.  2.  At 
AUV’s  top  or  bottom.  Either  upward  or  downward  flow  velocity  measurement  would 
be  in  the  wake,  and  the  AUV’s  launch  gear  or  the  storage  cart  would  jeopardize  the 
ADV  probe.  3.  At  the  horizontal-plane  flank  of  the  vehicle.  The  ADV  probe  would 
interfere  with  the  vehicle’s  recovery  hoop. 

Our  choice  is  to  mount  the  ADV  at  the  vehicle’s  largest  cross-section,  with  its 
probe  pointing  45°  from  the  vehicle’s  horizontal  central  plane,  as  shown  in  Figure  8.2. 
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Figure  8.2:  Cross-sectional  view  and  side  view  of  the  ADV’s  mounting  on  the  vehicle 
for  the.  1998  Labrador  Sea  Experiment  (in  the  lower  panel,  the  lower  half  of  the  AUV’s 
inner  fairing  is  placed  upside-down  for  ease  of  installation) . 

The  ADV’s  three  receiver  tips  reach  the  brink  of  the  hull  but  do  not  protrude  beyond 
it.  ■  This  satisfies  all  the  requirements  listed  above.  In  the  vehicle,  the  ADV  probe 
is  mounted  with  a  horizontal  plate  and  a  45°  slanted  bracket.  During  installation;  a 
laser  pointer  is  used  to  improve  precision  of  alignment. 

■  The  ADV  works  at  24  V  provided  by  the  AUV’s  batteries.  ADV  data  are  read 
by  the  vehicle’s  computer  through  a  serial  port,  and  saved  in  the  vehicle’s  central 
data  file.  This  guarantees  time  synchronization  with  other  instrument  data.  The 
bandwidth  of  the  AUV’s  data  bus  limits  the  ADV’s, data  output  rate  to  2,5  Hz, 
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8.2  Algorithm  of  Extracting  Earth- Referenced  Ver¬ 
tical  Flow  Velocity  from  AUV’s  Raw  Measure¬ 
ment 


The  AUV-borne  ADV  measures  flow  velocity  relative  to  the  moving  vehicle.  Hence  to 
obtain  the  Earth-referenced  flow  velocity,  i.e.,  the  true  flow  velocity,  we  must  subtract 
the  vehicle’s  own  velocity  from  the  raw  measurement.  Another  issue  of  concern  is  the 
vehicle  hull’s  influence  on  the  measurement.  As  shown  in  Figure  8.2,  the  ADV  probe’s 
sampling  volume  is  located  about  13  cm  from  the  vehicle’s  hull  surface  (the  vehicle  has 
a  length  of  2.2  m,  and  a  diameter  of  0.6  m  at  its  largest  vertical  cross-section  where 
the  ADV  is  mounted) .  This  distance  is  small  enough  to  necessitate  the  consideration 
of  the  hull’s  influence  on  flow  measurement.  In  this  section,  we  present  the  algorithm 
of  extracting  the  Earth-referenced  vertical  flow  velocity. 


u 


•►X(AUV) 


T 

Y(AUV) 

Figure  8.3:  Definition  of  velocity  vectors  in  the  AUV  coordinate  system  (plan  view). 

Figure  8.3  shows  the  plan  view  of  the  spheroid  that  serves  as  our  model  of  the 
AUV  [76].  We  choose  the  spheroid  to  be  of  the  same  aspect  ratio  as  the  Odyssey 
IIB  AUV’s  axisymmetric  fairing.  The  length  of  its  major  and  minor  axes  are  a 
and  0.293a,  respectively,  where  a=0.991  m.  This  represents  fairly  closely  the  shape 
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of  the  fore  half  of  the  hull,  which  poses  the  most  important  influence  on  the  flow 
measurement  as  the  vehicle  is  running  forward.  The  ADV’s  sampling  volume  is 
located  at  rm  =  [0  0.304a  0.304a]r.  At  a  normal  straight-ahead  flight  attitude  of 
the  AUV,  the  ADV’s  sampling  volume  is  far  outside  the  boundary  layer  or  any  wakes, 
in  accordance  with  the  mounting  criteria  laid  out  in  Section  8.1. 

The  AUV  coordinate  system  is  at  rest  in  the  Earth  coordinate  system,  but  its 
origin  and  orientation  is  coincident  with  the  vehicle’s  at  any  instant.  Velocities  U,  V, 
54,  vm,  and  um  listed  below  are  all  in  the  AUV  coordinate  system.  U  =  [ JJ\  U2  t/3]r 
is  the  current  velocity,  which  is  assumed  to  be  constant  over  the  length  of  the  vehicle. 
The  vehicle  translates  with  speed  V  =  [Vi  V2  V3]T,  and  rotates  with  angular  velocity 
54  =  [Oi  542  543]U  These  motions  are  also  referred  to  as  degree  of  freedom  1  through 
6,  respectively.  The  sampling  volume,  located  at  rm,  translates  with  velocity 


vm=V  +  fixrm  (8.2) 

The  vehicle’s  motion  relative  to  the  water  surrounding  it  imparts  a  disturbance 
at  rm.  The  flow  velocity  that  the  ADV  measures,  um,  therefore  has  a  difference 
from  U  —  vm.  According  to  the  potential  flow  theory,  the  disturbance  is  a  linear 
combination  of  the  components  of  V  —  U  and  54.  So  um  equals  U  —  vm  plus  a 
correction  term: 


um  =  U  -  vm  +  [A(V  -  U)  +  £54] 

=  (A-  /)(V  -  U)  +  (B  +  rmx)54 

where  Equation  (8.2)  is  incorporated.  A  and  B  are  two  square  matrices  describing 
the  AUV  hull’s  effect  on  flow  velocity  measurement:  A  for  translational  motion  and  B 
for  rotational  motion.  By  potential  flow  theory  calculations  conducted  by  Dr.  Knut 
Streitlien  [76],  the  two  matrices  are  found  to  be 
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0.0678  0 


0 


(8.4) 


A=  0  0.0563  0.4053 

0  0.4053  0.0563 

0  0.0687  -0.0687 
£=  0  0  0  a  (8.5) 

0  0  0 

Thus  for  instance,  if  the  vehicle  is  moving  forward  with  velocity  [1  0  0]T  rel¬ 
ative  to  the  ambient  current,  the  ADV’s  sampling  volume  will  see  a  velocity  of 
[—1.068  0  0]T  due  to  the  accelerated  flow  at  the  hull’s  maximum  diameter. 

The  Earth-referenced  flow  velocity  is  extracted  through  the  following  steps. 

1.  Transform  the  velocity  measurement  in  the  ADV  coordinate  system  into  the 
AUV  coordinate  system,  both  systems  shown  in  Figure  8.2. 

um  =  TadV->AUV  UADV  (8.6) 

where  uAdv  is  the  velocity  vector  originally  measured  in  the  ADV  coordinate 

system. 

Matrix 

0-10 

Tadv^auv  =  — sin(a )  0  — cos(a )  (8.7) 

cos(a)  0  —sin(a) 

transforms  from  the  ADV  coordinate  system  to  the  AUV  coordinate  system, 
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where  a  =  45°  is  the  mounting  angle  of  the  ADV  probe’s  stem  relative  to  the 
vehicle’s  horizontal  central  plane. 

2.  Compensate  the  AUV  hull’s  influence  and  subtract  the  velocity  induced  by  the 
vehicle’s  rotation,  by  applying  Equation  (8.3)  and  Equation  (8.6).  Then  the 
relative  flow  velocity  in  the  AUV  coordinate  system  is  obtained: 


U  -  V  =  (/  -  A)-1[um  —  (B  +  rmx)D] 

=  (I  —  A)  1{TADv->auv  uadv  —  (B  +  rmx)Q] 


(8.8) 


where  the  angular  velocity  (i.e.,  the  vehicle’s  yaw/pitch/roll  rate)  vector  Q,  is 
measured  by  the  AUV’s  KVH  Digital  Gyro  Compass  and  Digital  Gyro  Incli¬ 
nometer  [77]. 

3.  Recover  the  relative  flow  velocity  in  the  Earth  coordinate  system  using  the 
vehicle’s  heading,  pitch,  and  roll  measurements: 


U  Earth  V Earth  — 


U  l 
«2 


«3 


—  TAUV^Earth  (U-V)  (8.9) 

=  TAUV^-Earth  {I  ~  A)~1\Tadv^auv  Uadv  —  (B  +  rmx)f2] 


where 


TaUV^t  Earth  =  Th  Tp  Tr 


cos(eh ) 

-sin(9h) 

0 

cos(9p ) 

0 

sin(9p ) 

1 

0 

0 

sin(9h) 

cos(9h) 

0 

0 

1 

0 

0 

cos(9'r) 

-sin(9'r)^ 

0 

0 

1 

—sin(9p) 

0 

cos(9p ) 

0 

sin{9'r) 

cos(9’r ) 
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where  9h,  9P ,  and  9r  are  the  AUV’s  heading,  pitch,  and  corrected  roll  angles, 
respectively.  Transformation  from  the  Earth  coordinate  system  to  the  AUV 
coordinate  system  is  illustrated  by  Figure  8.4  which  contains  the  three  rotation 
angles.  Note  that  Equation  (8.10)  represents  a  transformation  in  the  oppo¬ 
site  direction:  from  the  AUV  frame  to  the  Earth  frame.  So  TAUV^Earth  = 
{TEarth-+Auv)~l ■  Since  the  three  rotational  matrices  are  all  orthogonal,  the  in¬ 
verse  of  T Earths auv  equals  the  product  of  transpositions  of  three  matrices  in 
the  reverse  order,  which  leads  to  Equation  (8.10). 


(X,Y,Z) :  E^rth  referarce  systan  - (x,y,z) :  AD V  referaxa  system 

Figure  8.4:  Transformation  from  the  Earth  coordinate  system  to  the  AUV  coordinate 
system. 

The  measured  roll  angle  9r  is  referenced  to  the  horizontal  plane,  rather  than  the 
rotational  angle  needed  in  matrix  Tr.  Trigonometric  derivation  [78]  gives  the 
relation  between  the  desired  rotational  angle  0r’  and  the  roll  measurement  9r  as 
expressed  by  Equation  (8.11).  When  the  pitch  angle  9P  is  small,  the  difference 
between  9r  and  9r ’  is  small  too. 


/  cos(29r)  +  cos(29v)  .  . 

9r  =  acos{- - — -4  v  p)  x  sign(9r ) 

2  cos(9p)  yJcos2{9p)  —  sin2(9r) 


(8.11) 


4.  Extract  the  Earth-referenced  flow  velocity,  i.e.,  the  true  flow  velocity.  In  the 
thesis,  we  are  concerned  about  the  vertical  component,  i.e.,  the  third  compo¬ 
nent  of  U Earth-  Let  us  denote  the  third  component  of  U EaTth  and  V Earth  as 
'tu  down -Earth  and  Vdown-Auv  i  respectively  (note  that  the  z-axis  points  downward, 
as  shown  in  Figure  8.4).  Then  according  to  Equation  (8.9),  we  have 
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Wdown-Earth  V down-AUV  —  ^.3 


(8.12) 


So  the  Earth-referenced  downward  flow  velocity  is 


^down^Earth  —  ^3  T  Vdown-AUV 

d  <8-13) 

=  u3  +  —  [AUV  depth  z(i)] 

(JLL 

where  m3  is  expressed  in  Equation  (8.9).  The  vehicle’s  own  vertical  velocity  is 
obtained  by  differentiating  its  depth  sensor  measurement. 

To  comply  with  the  commonly  adopted  upward  convention,  the  Earth-referenced 
upward  flow  velocity  is  written  as 


^ Earth  W doumJEarth 

i  (8-14) 

=  —  (u3  +  —[AUV  depth  z(f)]) 

8.3  Calibration  Experiment  at  the  David  Taylor 
Model  Basin 

To  ascertain  the  AUV  hull’s  influence  on  the  flow  velocity  measurement,  we  carried 
out  a  calibration  experiment  in  the  David  Taylor  Model  Basin  (DTMB)  in  August 
1997.  In  the  tow  tank,  a  complete  AUV  hull  equipped  inside  with  an  ADV  was  towed 
by  the  carriage  at  different  attack  angles  under  different  speed.  The  experiment 
required  a  large  tank  cross-section  to  minimize  influence  from  the  tank  boundaries, 
and  also  required  a  precise  speed  control.  The  cross-sectional  area  of  the  DTMB  tow 
tank  is  about  370  times  that  of  the  AUV.  The  tank  carriage  has  an  accuracy  of  about 
2  cm/s.  To  be  away  from  the  tank  boundaries  as  much  as  possible,  the  AUV  was 
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towed  along  the  central  line,  and  at  a  depth  of  about  2.8  m. 


8.3.1  Experiment  Design 

Let  us  present  the  experimental  set-up  in  relation  to  Equation  (8.3)  and  Figure  8.3. 
In  consideration  of  the  available  facilities,  we  do  not  attempt  to  generate  the  vehicle’s 
rotational  motion,  so  Q,  —  0  in  Equation  (8.3).  The  tank  water  is  still,  so  U  =  0. 
Then  Equation  (8.3)  is  simplified  to 


um  =  (A  -  7)V  (8.15) 

Hence  the  calibration  experiment  is  to  actually  validate  matrix  A  as  given  in  Equa¬ 
tion  (8.4).  This  matrix  depicts  the  AUV  hull’s  influence  on  the  flow  velocity  mea¬ 
surement,  which  is  induced  by  the  vehicle’s  translational  motion.  At  a  series  of  AUV 
speed,  ADV-measured  flow  velocities  are  to  be  compared  with  theoretical  predictions 
computed  by  Equation  (8.15)  using  matrix  A  in  Equation  (8.4).  To  test  out  various 
flow  orientations  relative  to  the  AUV,  we  need  to  enable  different  combinations  of  the 
vehicle’s  yaw  and  pitch  angles. 

The  experimental  structure  is  illustrated  in  Figure  8.5.  The  structure  is  composed 
of  three  parts:  a  rotating  bracket,  a  wedge,  and  a  hull  platform  (please  see  Appendix  B 
for  details).  The  upper  rectangular  plate  of  the  rotating  bracket  attaches  the  whole 
load  to  the  tow  tank  carriage.  Its  lower  circular  plate  connects  to  the  wedge  via  four 
bolts.  This  circular  plate  has  multiple  bolt  holes  to  permit  yaw  angles  of  0°,  5°,  10°, 
15°,  30°,  and  45°.  The  wedge  is  for  realizing  the  AUV’s  pitch/roll  angles  of  5°,  10°, 
15°,  30°.  The  hull  platform’s  upper  circular  plate  connects  to  the  wedge,  and  its 
lower  rectangular  plate  attaches  to  the  AUV’s  inner  fairing.  Its  45°  slanted  clamp 
(not  visible  in  Figure  8.5  from  this  perspective)  holds  the  ADV  probe. 

It  should  be  noted  that  in  this  calibration  experiment,  the  ADV  probe  points 
45°  upward  on  the  vehicle’s  port  side,  as  shown  in  Figure  8.5.  In  the  Labrador  Sea 
experiment,  the  ADV  probe  points  45°  downward  on  the  vehicle’s  starboard  side,  as 
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Figure  8.5:  Assembly  diagram  (upper)  for  the  DTMB  tow  tank  experiment;  and 
photo  (lower)  of  the  AUV  being  mounted  on  the  carriage  strut  using  a  5°  wedge,  with 
three  ADV  receiver  tips  visible.  •  - 

shown  in  Figure  8.2.  The  purpose  of  the  change  is  to  facilitate  recovery  of  the  AUV 
at  the  end  of  missions,  which  requires  contacts  at  the  upper  half  of  the  vehicle.  This 
difference  is  trivial  since  it  is  equivalent  to  rotating  the  vehicle  for  180°  about  its  along- 
ship  axis.  It  can  be  shown  that  this  reciprocal  move,  of  the  ADV  position  from  the  . 
port  side  to  the  starboard  side  causes  no  change  to  matrix  A,  and  only  sign  flippings 
in  Tadv^auv  (the  coordinate  transformation  matrix  given  in  Equation  (8.7)), 

In  the  1998  Labrador  Sea  Experiment,  the  AUV  had  a  “V”-shape  latch  at  its  noSe  ■ 
for  docking  to  an  underwater  station.  The  latch’s  length  is  about  30%  of  the  vehicle’s, 
and  its  largest  span  between  the  two  tips  equals  the.  vehicle’s  largest  diameter.  Its 
thickness  is  about  1  cm.  To  be  as  close  as  possible  to  field  operations,  we  also  added  a 
latch  to  the  vehicle  during  the  calibration  experiment.  It  is  found  that  the  latch  affects 
flow  velocity  measurement  only  at  “unfavorable”  attack  angles  when  the  vehicle  speed 
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reaches  2  knots.  This  will  be  shown  at  the  end  of  Subsection  8.3.2. 

As  shown  in  Figure  8.3,  the  AUV  coordinate  system  is  defined  such  that  its  x-axis 
points  forward,  y-axis  to  starboard,  and  z-axis  downward.  Accordingly,  yaw  is  the 
angle  between  the  AUV’s  along-ship  central  vertical  plane  and  that  of  the  tow  tank; 
pitch  is  the  angle  between  the  AUV’s  x-axis  and  the  horizontal  plane.  A  plus  sign 
of  yaw  means  that  the  AUV  steers  to  the  starboard  side,  while  for  pitch  it  means 
the  vehicle’s  nose  is  up.  Different  combinations  of  yaw  and  pitch  angles  as  shown 
in  Table  8.2  were  tested.  At  each  AUV  yaw/pitch,  the  carriage  ran  successively  at  1 
knot,  2  knots,  and  3  knots,  each  speed  lasting  for  about  40  s. 

Before  the  calibration  experiment,  we  distributed  neutrally  buoyant  hollow  glass 
spheres  [79]  in  the  tank  water.  The  tiny  spheres  (10~20  /tm)  acted  as  sound  scatterers, 
providing  strong  echoes  for  the  ADV’s  measurement.  There  are  two  parameters  for 
evaluating  the  ADV  data  quality:  the  correlation  coefficient  should  be  above  70% 
and  the  Signal  to  Noise  Ratio  (SNR)  should  be  above  10  dB,  to  achieve  the  1% 
velocity  precision  as  listed  in  Table  8.1.  In  the  calibration  experiment,  the  correlation 
coefficient  was  about  90%  and  the  SNR  was  about  15  dB. 

Two  installation  errors  are  calibrated  and  corrected  in  data  post-processing:  1. 
rotation  of  the  ADV  probe  in  the  clamp  (the  nominal  position  of  the  ADV’s  x- 
axis  should  be  in  the  vertical  plane),  2.  misalignment  between  the  hull  platform’s 
centerline  and  that  of  the  AUV.  The  first  error  is  found  to  be  2.7°  by  using  the  zero- 
yaw,  zero-pitch  calibration  run  without  the  AUV  hull.  The  second  error  is  found  to 
be  2.2°  by  using  a  zero- yaw,  zero-pitch  run  with  the  AUV  hull  on.  Note  that  the  2.2° 
error  angle  is  on  the  AUV’s  x-y  plane,  but  not  the  horizontal  plane  if  the  pitch  angle 
is  non-zero.  The  translation  formula  is 


,  tan(2.2  )  x 

yaw  error  =  tan  ( - - - - ) 

cos(pitch)' 


For  example,  at  the  largest  pitch  angle  -15°,  yaw  error=2.28°. 
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Table  8.2:  Tested  combinations  of  yaw  and  pitch  angles  (°) 


PitchYaw 

-15 

-5 

0 

5 

15 

-15 

X 

X  * 

-5 

X 

X 

X 

X 

WM 

0 

X 

X 

X  *  t 

X 

ma 

5 

Xf 

X 

X 

X:  ADV  Mounted  in  hull,  latch  attached, 
f:  ADV  Mounted  in  hull,  no  latch. 

*:  ADV  only. 


8.3.2  Comparison  of  Theoretical  and  Experimental  Results 

Figure  8.6  shows  the  comparison  of  experimental  results  (with  the  AUV’s  latch  on) 
and  the  theoretical  predictions  [76].  vx,  vy,  and  vz  are  in  the  ADV’s  coordinate 
system.  The  two  installation  errors  (given  at  the  end  of  Subsection  8.3.1)  have  been 
compensated.  Each  velocity  point  in  Figure  8.6  is  the  mean  value  of  600  data  points 
under  a  constant  carriage  speed,  which  is  further  normalized  by  the  carriage  speed 
for  display.  The  overall  measurement  noise  for  each  mean  velocity  is  reduced  to 
y  ( 2  cm/^s)26oQ1  cm/s)2  «  0.1  cm/s,  where  2  cm/s  is  the  carriage  speed  uncertainty 
and  1  cm/s  is  the  ADV’s  measurement  noise  (at  25  Hz  sampling  frequency  in  this 
experiment),  which  are  two  independent  errors.  In  the  above  calculation  it  is  assumed 
that  measurement  errors  for  different  data  points  are  uncorrelated. 

Relative  errors  between  the  experimental  results  and  the  theoretical  predictions 
are  shown  in  Table  8.3,  defined  as 


relative  error  = 


1 1 V experiment  F theory  2 

— } 

|  j  F experiment  |  1 2 


where  ||  •  ||2  denotes  the  Euclidean  norm.  Relative  errors  at  different  carriage  speed 
(1,  2,  3  knots)  are  averaged  to  give  the  tabulated  values. 

The  AUV  latch’s  influence  on  the  ADV’s  flow  measurement  is  shown  in  Figure  8.7. 
At  zero-yaw  and  zero-pitch,  the  influence  is  minimal:  measurements  with  and  without 
the  latch  are  very  close.  At  yaw=-5°  and  pitch=5°,  the  two  measurements  are  still 
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Pitch  =  -15°  '  Pitch  =  -5° 


Pitch  =  0° 


Pitch  =  5° 


Figure  8.6:  Comparison  of  the  calibration  experiment  results  and  the  theoretical 
predictions. 


very  close  at  carriage  speed  1  knot;  at  higher  speed  of  2  and  3  knots,  however,  the 
velocity  measurement  data  are  no  longer  valid  with  the  latch  on.  We  call  this  attitude 
“unfavorable”  for  the  ADV:  the  vehicle  headed  to  the  port  side  and  its  nose  was  up, 
so  the  ADV  (mounted  at  the  upper  port  side)  suffered  from  the  wake  generated  by 
the  protruding  latch,  although  the  effect  can  be  tolerated  at  low  speed  (1  knot). 

8.4  Error  Analysis  of  Processed  Vertical  Flow  Ve¬ 
locity 

The  total  estimation  error  of  the  extracted  Earth-referenced  vertical  flow  velocity 
w Earth  (expressed  in  Equation  (8.14))  has  two  components:  1.  noise,  2.  bias  due  to 
instrument  alignment  errors.  We  first  investigate  the  noise.  At  the  end  of  this  section, 
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Table  8.3:  Relative  errors  between  experimental  results  and  theoretical  predictions 

Pitch Yaw  1  -15°  -5°  0°  |  5°  |  15° 

■  ~T5°  4%  2% 


Pitch  Yaw 

-15° 

-5° 

0° 

5° 

-15° 

4% 

'  -5° 

;  4% 

3% 

3% 

3% 

0° 

3% 

2% 

2% 

3% 

-  5° 

3%* 

2% 

2% 

h  Only  at  i  knot  carriage  speed. 


yaw  =  0  ,  pitch  =  0 


yaw  =  -5  ,  pitch  =  5 


— 

6 

+ 

0 

-  Vx  (theory) 

-  Vy  (theory) 

Vz  (theory) 
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Figure  8.7:  Effects  of  AUV’s  latch  on  flow  measurement, 
we  will  discuss  the  bias. 

'  .The  estimation  noise  of  w Earth  results  from  three  sources  of  measurement  noise: 
i.  ADV,  2.  KVH  heading/pitch/roll  and  rate  sensor,  3.  AUV’s  depth  sensor.  We 
trace  the  evolution  of  noise  following  the  same  steps  as  in  Section  8.2: 

.1.  Consider  the  covariance  matrix  of  ADV’s  three-axis'  velocity  measurements  to 
be  Tjadv  —  °\dv  !■>  where. we  regard  noise  of  the  three  velocity  measurements 
as  uncorrelated.  Since  the  ADV’s  x  and  y-axis  have  higher  noise  than  z-axfs, 
we  uniformly  take  the  x-axis  measurement  noise  as  oadv •  Then  through  the 
orthogonal  transformation  in  Equation  (8.6),  the  covariance  matrix  of  um  is 


Coi^Um]  —  TADV-^AUV  &ADV  ^  ^ADV^AUV 

_  2  T 

~  °ADV  1 

where  we  note  that  Tadv-+auv  is  an  orthogonal  matrix. 

2.  By  Equation  (8.8), 


(8.16) 


Cov[ U  -  V]  =  (/  -  A)~'Cov\ um  -  (B  +  rmx)Sl]((/  -  /t)-‘)T 

(8.17) 

Noting  that  the  noise  of  ADV  measurement  um  is  independent  of  that  of  KVH 
heading/pitch/roll  rate  measurement  Q,  we  have 


Cov[ um  -  (B  +  rmx)f2]  =  Cou[um]  +  Cov[(B  +  rmx)Q] 

=  Cov{  um]  +  Cov[(B  +  R)tt} 


(8.18) 


where  B  is  expressed  in  Equation  (8.5)  and 


0  -rz  ry 

rz  0  —rx 

-ry  rx  0 


(8.19) 


whose  elements  are  the  three  components  of  rm  —  the  location  vector  of  the 
ADV’s  sampling  volume,  as  given  at  the  beginning  of  Section  8.2. 

Consider  the  covariance  matrix  of  the  KVH’s  rate  measurements  to  be  = 
ctq  I.  The  pitch/roll  rate  measurement  noise  is  slightly  larger  than  that  of  the 
heading  rate,  so  we  uniformly  take  the  larger  one  as  gq,.  Then  by  incorporat- 
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ing  Equation  (8.16),  Equation  (8.18)  becomes 


Cov[nm  -  (B  +  rmx)fi]  =  a\DV  I  +  a^(B  +  R)(B  +  R)T 

(8.20) 


Thus  Equation  (8.17)  becomes 


Cov[ U  -  V]  =  (/  -  A)-'[a\DV  I  +  a2n(B  +  R)(B  +  -  7l)->)T 

(8.21) 


3.  By  Equation  (8.9),  we  have 


Cov[U  Earth  ^  Earth] 

=  Cov{[ui  u2  u3]t}  (8.22) 

=  R AUV ^tEarth  Cov\ U  —  V]  T^uv_^Earth 

where  TAuv^Earth  —  Th  Tp  Tr  as  in  Equation  (8.10). 

4.  By  Equation  (8.14),  we  have 


aWBarth  °M3  °VUV_AVV  (8.23) 

where  vup_auv  =  —Vdown^Auv  (in  Equation  (8.13)),  computed  by  differentiat¬ 
ing  the  vehicle’s  depth  reading.  Note  that  the  noise  of  the  calculated  u3  is 
uncorrelated  with  the  AUV’s  depth  measurement  noise. 

By  definition,  <j„3  is  element  (3,3)  of  Cov[V Earth  —  V Earth]  in  Equation  (8.22), 
so 
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G WEarth  {element  (3,  3)  of  Cov\\J  Earth  "V Earth]}  "f"  & up-AUV 

(8.24) 

In  level  runs  as  we  will  analyze  for  AUV  Mission  B9804107,  the  vehicle’s  pitch 
and  roll  angles  are  close  to  zero.  Hence  when  calculating  Cov[UEarth  ~  VEQrih] 
using  Equation  (8.22),  we  can  approximately  deem  Tp  =  Tr  —  I  (please  refer 
to  Equation  (8.10)).  Then  TAUV_+Earth  ~  Th.  Since  the  third  row  of  Th  is  [0  0  1], 
TAuv^Earth  Cov[ U  -  V]  TAuv^Earth  preserves  element  (3,3)  of  Cov{ U  -  V], 
Thus  by  incorporating  Equation  (8.21)  and  Equation  (8.22),  Equation  (8.24) 
becomes 


w  Earth 

~  {element  (3, 3)  of  Cov (U  -  V]}  +  a2Vup_AUV 

=  {element  (3, 3)  of  (/  -  A)~x[o2adv  I  +  o2u{B  +  R){B  +  R)T]{(I  -  A)~l)T}  + 

Now  let  us  quantify  the  three  error  sources.  In  the  following,  the  values  of  error 
No.  2  (of  KVH  measurement)  and  error  No.  3  (of  depth  measurement)  do  not  vary 
with  AUV  missions.  The  value  of  error  No.  1  (of  ADV  measurement)  is  specifically  for 
AUV  Mission  B9804107,  which  will  be  analyzed  in  detail  in  Chapter  9.  The  ADV’s 
measurement  noise  depends  on  operation  conditions  in  different  regions  of  water. 

•  For  AUV  Mission  B9804107,  there  appears  a  noise  floor  of  5  cm/s  in  the  pe- 

riodograms  of  ADV’s  x  and  y  velocity  measurements  (z  velocity  is  less  noisy). 
The  ADV’s  data  output  rate  is  2.5  Hz.  Then  50-s  smoothing  is  equivalent  to 
125-point  averaging.  Assuming  white  noise,  oAdv  =  «  0.45  cm/s  after 

smoothing. 

•  The  KVH  pitch/roll  sensor’s  repeatability  is  0.25°  (heading’s  is  smaller)  [77]. 
The  KVH  sensor’s  data  rate  is  5  Hz.  So  by  differentiation,  noise  of  rate  is 
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0.25°  x  vr/ 180' 
0.2  s 


=  0.022  rad/s.  50-s  smoothing  is  equivalent  to  250-point  averaging. 
Thus  op  =  °'02jTTrf/,s  1.4  x  1CV3  rad/s  after  smoothing. 

•  The  AUV’s  depth  sensor  (Model  8B-4000)  [80]  has  a  precision  of  0.01  m.  The 
depth  sensor’s  data  rate  is  about  5  Hz.  So  by  differentiation,  the  noise  of 
the  vehicle’s  vertical  velocity  estimate  is  °0°*  ™  =  5  cm/s.  50-s  smoothing  is 
equivalent  to  250-point  averaging.  Thus  a „uMirv  =  ~  0.32  cm/s  after 

smoothing. 

Incorporating  the  above  ctadv ,  op,  avup_AUV,  as  well  as  matrix  A  (Equation  (8.4)), 
B  (Equation  (8.5)),  and  R  (Equation  (8.19))  into  Equation  (8.25),  we  get  crWEarih  ~ 
0.7  cm/s.  These  results  are  summarized  in  Table  8.4. 


Table  8.4:  Measurement/estimation  noise 


Measurement/ 

estimation 

Flow 

velocity 

Heading/ attitude 
rate 

AUV  vertical 
velocity 

wEarth 

Sensor 

ADVOcean 

5  MHz 

KVH 

Paroscientific 

8B-4000 

Noise  after 
50-s  smoothing 

0.45  cm/s* 

1.4  x  10“3  rad/s 

0.32  cm/s 

0.7  cm/s* 

V  For  AUV  M 

ission  B9804107  in  the  1998  Labrador  Sea  Experiment. 

Besides  measurement  noise,  bias  is  another  type  of  error.  In  the  estimation  of 
w Earth >  installation  errors  of  the  ADV  probe  and  the  KVH  box  are  of  concern.  In 
processing  the  data  for  AUV  Mission  B9804107,  two  errors  are  corrected  as  shown 
in  Figure  8.8:  2°  rotation  of  the  ADV  probe  (face  view)  and  1°  pitch-up  of  the  KVH 
box.  Although  we  used  a  laser  pointer  to  improve  alignment  precision,  accurate 
installation  of  the  ADV  probe  was  still  difficult.  A  photograph  taken  by  facing  the 
probe  suggests  a  clockwise  error  as  shown  in  Figure  8.8.  The  KVH  box  is  mounted 
in  the  front  sphere  of  the  AUV.  Precise  alignment  is  also  very  challenging,  both  for 
the  KVH  box  inside  the  sphere  and  the  glass  sphere  itself.  The  values  of  errors  are 
hard  to  ascertain.  The  criterion  we  take  in  data  processing  is:  after  error  corrections, 
the  mean  value  of  convective  WEarth  during  an  AUV  level  run  should  vanish.  This 
criterion  is  based  on  mass  preservation  and  also  the  output  data  of  the  convection 
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model  (presented  in  Subsection  5.1.4).  For  AUV  Mission  B9804107,  after  correcting 
the  above  two  installation  errors,  wEarth  at  the  20-m  depth  is  about  -0.2  cm/s.  w Earth 
at  the  250-m  depth  is  about  -0.9  cm/s,  but  this  larger  magnitude  is  largely  due  to 
a  negative  mean  during  the  third  and  fourth  legs  (in  Section  9.6,  we  will  discuss 
possible  reasons).  With  improvement  of  installation  accuracy,  the  problem  of  bias 
will  be  reduced. 

1  deg. 

AUV  bow  -  . 


2  deg. 


Figure  8.8:  Installation  errors  that  are  corrected  in  estimating  wEarth  in  AUV  Mission 
B9804107. 
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Chapter  9 


Labrador  Sea  Experiment 

9.1  Background 


Cruise  Track,  Kn156 


Figure  9.1:  Ship  track  of  R/V  Knorr  during  the  1998  Labrador  Sea  Experiment 
(courtesy  of  Dr.  Knut  Streitlien) . 


The  Labrador  Sea  lies  between  northern  Canada  and  Greenland.  It  is  one  of  the 
few  locations  in  the  world  where  open  ocean  convection  occurs  [45],  [46].  During 
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the  winter,  the  sea  surface  is  subjected  to  intense  heat  flux  to  the  atmosphere.  The 
resulting  buoyancy  loss  causes  the  surface  water  to  sink  to  large  depths,  initiating 
ocean  convection.  Detailed  discussions  are  found  in  Subsection  5.1.1. 

During  January/February  1998,  researchers  from  MIT,  the  Woods  Hole  Oceano¬ 
graphic  Institution,  and  the  University  of  Washington,  led  an  expedition  to  the 
Labrador  Sea  to  study  ocean  convection.  The  Research  Vessel  (R/V)  Knorr  was 
employed  in  this  expedition.  The  map  of  the  Labrador  Sea  area  as  well  as  the  ship 
track  is  shown  in  Figure  9.1  (The  cruise  number  was  KN156).  AUVs  and  other 
oceanographic  platforms  were  deployed  in  this  experiment.  The  current  flow  velocity 
was  measured  by  an  Acoustic  Doppler  Velocimeter  (ADV)  installed  in  an  Odyssey  IIB 
AUV ,  as  presented  in  Section  8.1.  Vertical  flow  velocity  is  a  crucial  signature  of  ocean 
convection  [47],  [16].  It  is  the  quantity  we  use  in  the  thesis  for  distinguishing  ocean 
convection  from  internal  wave.  In  the  ending  section  of  this  chapter,  the  AUV-based 
classifier  will  be  tested  by  the  Labrador  Sea  data. 


9.2  Observed  Conditions  for  Ocean  Convection 

During  the  January/February  1998  Labrador  Sea  Experiment,  the  meteorological 
measurement  recorded  intense  heat  loss  from  the  ocean  to  the  atmosphere.  Hydro- 
graphic  measurements  by  the  AUV  and  on-deck  CTD  casts  showed  mixed  water  layer 
down  to  several  hundred  meters.  Intense  surface  heat  loss  plus  vertically  mixed  water 
layer,  provides  very  good  conditions  for  ocean  convection.  In  Section  9.3,  we  will 
introduce  direct  observations  of  convection  using  Lagrangian  floats  of  the  University 
of  Washington. 

9.2.1  Meteorological  Condition 

Meteorological  data  were  recorded  by  an  Improved  Meteorological  (IMET)  system  [81] 
on  board  R/V  Knorr.  Prof.  Peter  Guest  of  the  Naval  Postgraduate  School  calculated 
the  ocean  surface  heat  flux  based  on  the  measurements.  The  time  series  of  total 
surface  heat  flux  is  shown  in  the  third  panel  of  Figure  9.2.  Note  that  by  the  convention 
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Figure  9.2:  Total  surface  heat  flux  time  series  (the  third  panel)  during  the  Jan¬ 
uary/February  1998  Labrador  Sea  Experiment  (courtesy  of  Prof.  Peter  Guest).. 

used  in  this  calculation,  negative  heat  flux  means  that  the  sea  water  loses  heat  to  the 
air.  ■  • 

‘  .  The  average  heat  flux  value  was  about  220  W/m2.  Around  January  24  and  Febru¬ 
ary  16,  the  heat  flux  value  was  as  large  as  600  W/m2  and  800  W/m2,  respectively. 
Meteorological  parameters  during  AUV  Mission  B98041Q7  (which  will  be  used  to  test 
the  classifier)  are  summarized  in  Table  9.1. 

Let  us  make  comparisons  with  previous  open  ocean' convection  experiments.  In 
the  Greenland  Sea  Experiment  [47]  during  the  winter  of  1988/1989,  the  heat  flux  fluc: 
tuated  between  100  W/m2  and  500  W/m2,  with  an  average  value  of  about  250  W/m2. 
Ocean  convection  was  observed  during  that  experiment,  using  moored  Acoustic  Doppler 

163 


Table  9.1:  Meteorological  parameters  (using  Prof.  Peter  Guest’s  calculation  results) 
during  AUV  Mission  B9804107 


Sea  surface  temperature 

Air  temperature 

Dewpoint 

Wind  speed 

Total  heat  flux* 

3.1°C 

-2.2°C 

-8.3°C 

10.3  m/s 

-300  W/m2 

*:  sign  means  that  the  sea  water  loses  heat  to  the  air. 


Current  Profilers  (ADCPs)  (with  30-minute  data  interval).  In  an  earlier  Labrador 
Sea  experiment  [46]  during  the  winter  of  1994/1995,  the  average  heat  flux  was  about 
300  W/m2.  Using  a  moored  ADCP  (with  20-minute  data  interval)  and  Profiling  Au¬ 
tonomous  LAgrangian  Circulation  Explorer  (PALACE)  floats,  ocean  convection  was 
observed.  The  sea  surface  heat  flux  value  in  our  Labrador  Sea  Experiment  is  close 
to  that  of  the  two  previous  experiments.  We  therefore  have  reason  to  expect  ocean 
convection’s  occurring  during  Cruise  KN156. 

9.2.2  Hydrographic  Condition  Measured  by  the  AUV  and 
On-Deck  CTD  Casts 

Besides  surface  heat  loss,  a  vertically  mixed  water  column  is  another  key  condition 
for  ocean  convection.  Across  the  Labrador  Sea  basin  (about  600  km),  a  mixed  water 
layer  of  depth  270  m  ~  500  m  was  observed  by  a  series  of  Conductivity-Temperature- 
Depth  (CTD)  casts  from  the  ship  deck.  During  two  different  AUV  missions  four  days 
apart,  mixed  water  layers  were  also  clearly  recorded.  During  Mission  B9804107,  the 
mixed  layer  was  down  to  350  m.  During  an  earlier  Mission  B9803703,  the  mixed  layer 
existed  over  a  horizontal  distance  of  about  8  km,  depth  varying  from  280  m  to  180  m. 

1.  On-deck  CTD  casts  across  the  Labrador  Sea  basin. 

From  the  ship  deck  of  Knorr,  Prof.  Eric  D’Asaro  and  his  group  conducted  CTD 
casts  at  20  stations  across  the  Labrador  Sea  basin,  as  shown  in  Figure  9.3.  This 
cross-section  is  known  as  “AR7W”  in  the  oceanographic  community  [49],  [45], 
covering  a  distance  of  over  600  km.  To  the  east  of  CTD  Station  No.  17  is  the 
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Figure  9.3:  Track  of  on-deck  CTD  cast  stations  across  the  Labrador  Sea  basin  (Section 
AR7W)  (Courtesy  of  Prof.  Eric  D’Asaro  and  Ms,  Elizabeth  Steffen). 

site  of  the  former  Ocean  Weather  Station  Bravo  (56°30’N,  51°00’W  [82]).- 

The  vertical  cross-section  of  potential  temperature  obtained  from  these  CTD 
casts  is  shown  in  Figure  9.4.  Corresponding  to  Figure  9.1,  the  left  boundary 
in  Figure  9.4  is  northern  Canada  and  the  right  boundary  is  Greenland.  Between 
CTD  Station  No.  16  and  No.  21,  the  water  was  well  mixed  down  to  about 
270  m  (measurements  made  on  January  31  ~  February  5  [83]).  Out  of  all  CTD 
stations  along  Section  AR7W,  these  six  were  closest  to  sites  of  AUV  missions  (to 
be  presented) .  The  observed  270-m  mixed  layer  depth  is  close  to  that,  measured 


Figure  9.4:  Potential  temperature  across  the  Labrador  Sea  basin  (Section  AR7W) 
(Courtesy  of  Prof.  Eric  D’Asaro  and  Ms.  Elizabeth  Steffen). 

by  the  AUV  during  Mission  B9803703,  to  be  given  below. 

Between  CTD  Station  No.  24  and  No.  27,-  however,  the  mixed  layer  was 
deeper:  about  500  m  (measurements  made  on  February  6  ~  February  8  [83]). 
Thus  we  see  a  large  variation  of  mixed  layer  depth  over  the  basin  scale  (several 
'  .  hundred  km).  The  impact  of  this  variation  on  AUV-based  classifier  design  will 
be  reviewed  at  the  end  of  this  section.  ■. 

2.  ' Mixed  layer  measured  by  the  AUV  during  Mission  B9803703. 
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In  contrast  to  the  above  basin-scale  hydrographic  sampling,  AUV  carried  out 
fine-scale  surveys  in  a  local  area.  AUV  Mission  B9803703  took  place  at  1:23 
~  3:30  on  February  6,  1998  (GMT).  The  mission  launch  location  was  about 
56°42’N,  52°46’W  —  where  the  AOSN  (Autonomous  Oceanographic  Sampling 
Network  [84])  mooring  was  anchored.  This  was  a  long-range  yo-yo  run  between 
190-m  and  400-m  depths,  covering  a  horizontal  distance  of  about  8  km. 

The  CTD  sensor  suite  on  the  vehicle  made  measurements  of  temperature,  salin¬ 
ity,  and  depth.  We  then  derive  potential  temperature  and  potential  density 
from  these  in-situ  measurements,  so  as  to  remove  water  pressure’s  effect.  Poten¬ 
tial  temperature,  salinity,  potential  density,  as  well  as  vehicle  depth  are  shown 
in  Figure  9.5. 


Figure  9.5:  Potential  temperature,  salinity,  potential  density,  and  AUV  depth  during 
Mission  B9803703. 

To  facilitate  inspection,  we  plot  potential  density  profiles  for  successive  de¬ 
scent/ascent  legs  during  this  yo-yo  mission.  In  regions  cruised  in  the  early 
stage,  water  was  mixed  down  to  about  280  m.  The  mixed  layer  became  shal- 
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lower  later  on,  up  to  about  180  m  where  the  vehicle  ended  the  mission.  As 
observed  by  the  vehicle  within  the  mission  range,  the  mixed  layer  extended  for 
8  km  with  a  thickness  variability  of  100  m. 
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Figure  9.6:  Potential  density  measured  on  successive  descent/ascent  legs  during  yo-yo 
Mission  B9803703  (density  anomaly,  i.e.,  potential  density  - 1000  kg/m3,  is  displayed). 

3.  Mixed  layer  measured  by  the  AUV  during  Mission  B9804107. 

AUV  Mission  B9804107  took  place  at  4:46  ~  7:00  on  February  10,  1998  (GMT), 
four  days  later  than  Mission  B9803703.  The  mission  location  was  about  the 
same  as  that  of  Mission  B9803703.  The  AUV  behaviors  in  this  mission  are 
illustrated  in  Figure  9.7.  The  vehicle  first  spiraled  down  to  426-m  depth,  then  it 
spiraled  up  to  250-m  depth.  At  this  depth  plane,  the  vehicle  made  a  “diamond” 
run,  i.e,  closed  a  4-leg  loop  with  90°  turns,  each  leg  lasting  for  720  s.  After  that, 
it  spiraled  up  to  20-m  depth,  making  an  identical  “diamond”  run.  At  the  end, 
the  vehicle  ascended  to  the  sea  surface. 

Profiles  of  potential  temperature,  salinity,  and  potential  density  are  shown 
in  Figure  9.8.  These  profiles  demonstrate  that  water  was  well  mixed  down 
to  350  meters.  Water  was  stratified  below  this  mixed  layer. 
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Figure  9.7:  AUV  behavior  sequence  in  Mission  B9804107. 


We  have  not  only  observed  significant  surface  cooling  and  mixed  water  layers,  but 
also  noticed  that  the  associated  parameters  (i.e.,  heat  flux  and  mixed  layer  depth) 
experience  variations  over  time  and  space.  This  necessitates  that  the  AUV-based  clas¬ 
sifier  should  be  robust  to  parameter  uncertainty.  In  the  robustness  analysis  presented 
in  Chapter  7,  we  treat  the  variation  of  these  key  parameters  as  the  “global  uncer¬ 
tainty”  .  Its  effect  is  to  increase  the  total  variance  of  PSD  estimate,  and  accordingly 
suppress  the  amplitude  of  the  feature  transformation  vector.  The  classifier’s  perfor¬ 
mance  (probability  of  detection  versus  probability  of  false  alarm)  is  thus  lowered,  but 
with  a  gain  of  robustness  to  model  mismatch. 

9.3  Convection  Shown  by  Lagrangian  Floats  of  the 
University  of  Washington 

Prof.  Eric  D’Asaro  of  the  University  of  Washington  developed  a  type  of  Lagrangian 
float  [10]  that  follows  water  motion  through  density  matching  and  a  high  drag.  Its 
horizontal  motion  is  tracked  acoustically  and  its  vertical  motion  is  recorded  by  pres¬ 
sure  measurement. 

In  the  1998  Labrador  Sea  Experiment,  seven  Lagrangian  floats  were  also  deployed 
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Potential  temperature  ( °  C)  Salinity  (psu)  Potential  density  (kg/m3) 

Figure  9.8:  Profiles  of  potential  temperature,  salinity,  and  potential  density  during 
AUV  Mission  B9804107. 

to  study  convection.  They  were  launched  on  January  26  at  56°30’N,  51°00’W  [85],  to 
the  northwest  of  AUV  mission  sites  (Mission  B9803703  and  B9804107).  The  floats’ 
depth-time  records  in  February  1998  are  shown  in  the  lower  panel  of  Figure  9.9.  A 
convectively  mixed  layer  extending  from  the  surface  to  about  600-m  depth  is  evident. 
The  floats  cycled  across  this  layer,  following  water’s  convective  motion.  The  time 
of  records  shown  in  Figure  9.9  was  two  weeks  (or  more)  later  than  the  two  AUV 
missions:  B9803703  (2/6/1998)  and  B9804107  (2/10/1998).  As  the  winter  progresses, 
the  convective  layer  continues  to  deepen  due  to  sustained  surface  cooling.  This  should 
explain  why  the  mixed  layer  depth  recorded  by  the  floats  is  larger  than  that  measured 
by  the  AUV.  The  Lagrangian  floats’  records  confirm  not  only  the  existence  of  mixed 
layers,  but  also  the  occurrence  of  convection. 

Furthermore,  the  root-mean-square  (rms)  vertical  flow  velocity  is  found  to  be 
2  ~  3  cm/s  based  on  the  float  data  [83].  In  the  AUV  data  processing  results  to  be 
shown  in  Figure  9.11  and  Figure  9.12,  the  counterpart  is  2  cm/s  at  the  250-m  depth 
and  3.5  cm/s  at  the  20-m  depth  (both  of  Mission  B9804107).  Measurements  by  these 
two  independent  platforms  are  consistent. 
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Figure  9.9:  Depth  history  of  Lagrangian  floats  in  the  Labrador  Sea,  in  1997  (top)  and 
1998  (bottom).  The  abscissa  is  day  of  February.  (Courtesy  of  Prof.  Eric  D’Asaro.) 

9.4  Data  Processing  Results  of  AUV-Measured  Ver 
tical  Flow  Velocity 

AUV  Mission  B9804107  has  been  described  in  Section  9.2.  The  mission  comprises 
two  long-duration  runs  at  20-m  and  250-m  depths,  as  shown  in  Figure  9.7.  On  the 
vehicle,  the  ADV  probe  measured  three-dimensional  flow  velocity,  while  the  CTD 
sensor  suite  measured  other  water  properties.  These  measurements  are  shown  in  Figr 
ure  9.10,  where  potential  temperature  and  potential  density  are  deduced  from  in-situ 
measurements.  • .  . 

In  the  Labrador  Sea,  the  AUV-borne  ADV  experienced  low  correlation  coefficient 


3-D  vel.  (m/s)  meas.  by  ADV  (after  50  s  smoothing)  along  with  CTD  and  potential  density  in  Mission  B9804107 
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Figure  9.10:  During  AUV  Mission  B9804107,  flow  velocity  was  measured  by  an  Acous¬ 
tic  Doppler  Velocimeter  (ADV).  Conductivity,  temperature,  and  depth  were  measured 
by  on-board  CTD  sensors. 


and  low  Signal  to  Noise  Ratio  (SNR).  This  was  probably  due  to  the  low  density  of 
floating  particles  in  the  water.  Measures  of  electrical/magnetic  noise  reduction,  such 
as  twisting  the  ADV’s  signal  wire  pairs,  will  also  be  necessary  for  improving  the  SNR. 

By  applying  a  50-s  smoothing  to  the  Labrador  Sea  data,  the  final  estimation  noise 
of  the  extracted  Earth-referenced  vertical  flow  velocity  w Earth  is  reduced  to  about 
0.7  cm/s,  as  shown  in  Table  8.4. 

In  AUV  Mission  B9804107,  an  acoustic  communication  device  (MOdulator/DEModulator, 
i.e.,  MODEM)  periodically  transmitted  high-power  signals  with  a  duty  cycle  of  3.5  s 
every  28  s.  In  the  3.5-s  transmission  duration,  the  ADV  data  were  found  to  be  cor¬ 
rupted  by  the  MODEM’S  interference,  probably  through  the  power  supply  system1. 

1In  a  later  AUV  experiment  in  1999  using  a  different  flow  measurement  sonar  —  an  Acoustic 
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Figure  9.11:  The  Earth- referenced  vertical  flow  velocity  w Earth  (the  first  panel)  at 
the  250-m  depth  of  AUV  Mission  B9804107.  In  the  second  panel  is  the  AUV’s  own 
vertical  velocity.  The  AUV’s  roll  (the  third  panel)  shows  when  the  vehicle  made  90° 
turns. 

We  apply  the  algorithm  in  Section  8.2  to  process  the  ADV  measurements  made 
at  250-m  and  20-m  depths  during  Mission  B9804107.  The  raw  measurements  are 
plotted  in  the  first  three  panels  in  Figure  9.10.  The  data  processing  results  are  shown 
in  Figure  9.11  and  Figure  9.12,  where  the  vehicle’s  own  vertical  velocity  (the  second 
panel)  has  been  removed  for  producing  the  Earth-referenced  vertical  flow  velocity 
W Earth  (the  first  panel). 


Doppler  Current  Profiler  (ADCP),  severe  interference  from  the  MODEM  transmission  was  again 
noted.  This  problem  needs  to  be  solved. 

2The  exclusion  ratio  in  data  processing  is  16%  due  to  the  leading  and  trailing  edges  of  the 
MODEM  pulses.  The  ADV-measured  along-ship  flow  velocity  (i.e.,  the  vehicle’s  cruise  speed)  is 
accordingly  reduced  after  the  50-s  smoothing,  but  the  exclusion  has  little  impact  on  the  transverse 
flow  velocity  (including  the  vertical  flow  velocity)  which  is  basically  zero-mean. 
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Figure  9.12:  The  Earth-referenced  vertical  flow  velocity  w Earth  (the  first  panel)  at 
the  20-m  depth  of  AUV  Mission  B9804107.  In  the  second  panel  is  the  AUV’s  own 
vertical  velocity.  The  AUV’s  roll  (the  third  panel)  shows  when  the  vehicle  made  90° 
turns. 

In  Figure  9.13,  we  overlap  WEarth  at  the  two  depths:  250-rn  and  20-m.  A  shift  of 
3550  s  is  added  to  the  time  index  for  the  250-m  depth,  during  which  the  vehicle  com¬ 
pleted  the  diamond  loop  at  the  250-m  depth  and  ascended  to  the  20-m  depth  plane 
(please  refer  to  Figure  9.7).  We  note  that  for  about  600  s  in  Figure  9.13,  consistency 
between  the  two  curves  is  pronounced.  According  to  the  MIT  Ocean  Convection 
Model  data  presented  in  Subsection  5.1.4  and  Subsection  6.2.1,  over  a  depth  differ7 
ence  of  230  m  and  a  time  difference  of  3550  s,  vertical  velocities  at  the  same  horizontal 
position  should  show  strong  similarity.  The  reason  is  that  convective  cells  are  verti¬ 
cally  aligned,  and  the  field  varies  very  slowly  as  convection  approaches  a  stationary 
state  (as  discussed  in  Subsection  6.2.1).  As  will  be  pointed  out  in  Section  9.6,  the 
convection  model  is  inadequate  for  depicting  the  20-m  depth  shallow  water  as  it  does 
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not  include  effects  of  the  atmospheric  forcing.  So  the  model-based  “strong  similar¬ 
ity”  between  depths  may  not  apply.  Nevertheless,  it  still  offers  a  means  for  checking 
measurements.  The  consistency  of  wEarth  between  the  two  depths  as  shown  in  Fig¬ 
ure  9.13,  although  only  lasting  for  a  short  duration,  adds  to  our  confidence  in  the 
data  processing  results.  , 


Overlapped  Earth-referenced  vertical  flow  velocities  at  250-m  and  20-m  depths 


Figure  9.13:  Comparison  of  wEarth  of  250-m  and  20-m  depths,  with  a  time  shift  of 
3550  s  corresponding  to  AUV’s  cruise  time  (please  see  Figure  9.7). 

9.5  Tests  of  AUV-Based  Classifier  Using  Labrador 
Sea  AUV  Data 

As  shown  in  Table  5.1  in  Subsection  5.1.3,  parameters  of  the  MIT  Convection  Model 
are  set  to  measured  values  (heat  flux,  mixed-layer  depth,  and  Coriolis  frequency  cor¬ 
responding  to  the  experiment  site)  during  AUV  Mission  B9804107.  The  model  output 
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of  vertical  flow  velocity  (one  snapshot  shown  in  Figure  5.2)  are  then  used  to  compute 
the  temporal-spatial  PSD  (i.e.,  the  rj-v  spectrum),  as  presented  in  Subsection  6.2.1. 
From  the  r]-v  spectrum,  we  obtain  the  corresponding  mingled  spectra  at  a  series  of 
AUV  speed,  as  shown  in  Figure  6.5.  These  mingled  spectra  constitute  the  theoreti¬ 
cal  template  for  convection.  The  AUV-based  classifier  is  built  upon  this  convection 
template  along  with  the  internal  wave  template,  as  presented  in  Section  6.3. 

In  Section  6.5,  we  have  tested  the  classifier  using  model-based  synthesized  data.  In 
this  section,  we  test  the  classifier  using  the  Labrador  Sea  data  acquired  by  the  AUV- 
borne  ADV.  The  input  to  the  classifier  is  the  Earth-referenced  vertical  flow  velocity 
(denoted  WEarth)  extracted  from  raw  data  during  AUV  Mission  B9804107,  using  the 
algorithm  in  Section  8.2  and  shown  in  Figure  9.11  (250-m  depth)  and  Figure  9.12 
(20-m  depth). 

9.5.1  At  250-m  depth  of  AUV  Mission  B9804107 

The  convection  model  parameters  use  measurements  in  this  AUV  mission.  Further¬ 
more,  model  computations  are  carried  out  at  the  depth  of  250  m.  We  therefore  expect 
to  see  that  the  model-based  classifier  recognizes  the  250-m  depth  AUV  data  as  con¬ 
vection.  The  PSD  estimate  of  wEarth  on  the  first  and  second  legs3  at  the  250-m  depth 
is  shown  in  Figure  9.14,  using  five-point  frequency-domain  smoothing  to  reduce  the 
variance  of  estimate.  We  note  a  spectral  peak  at  0.007  Hz.  The  AUV  speed  during 
this  mission  is  about  1  m/s.  As  shown  in  the  first  panel  of  Figure  6.5,  the  peak  fre¬ 
quency  of  the  mingled  spectrum  template  (based  on  the  convection  model)  for  AUV 
speed  1  m/s  lies  at  about  0.005  Hz.  These  two  frequencies  are  close. 

Now  let  us  feed  the  250-m  depth  wEarth  data  into  the  classifier.  The  PSD  estimate 
of  WEarth  is  shown  in  Figure  9.15,  along  with  the  two  PSD  templates  for  convection 
and  internal  wave  as  we  have  seen  in  the  upper  panel  of  Figure  6.6.  In  Figure  9.15, 
we  do  not  conduct  frequency-domain  smoothing  to  ensure  that  individual  frequency 
points  provide  uncorrelated  PSD  estimates  as  required  by  the  classifier’s  formulation 

3The  third  and  fourth  legs  will  be  discussed  in  Section  9.6. 
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PSD  of  the  1st  and  2nd  legs,  250-m  depth.  After  5-point  freq.  smoothing  with  A  f=7  x  ICf4  Hz. 


Figure  9.14:  PSD  estimate  of  WEarth  at  the  250-m  depth  of  AUV  Mission  B9804107. 
The  1-a  error  band  is  shown,  using  five-point  frequency-domain  smoothing. 

(detailed  discussions  are  found  at  the  end  of  Section  6.3). 

Processed  in  the  same  way  as  for  simulations  in  Section  6.5,  the  experimental 
test  result  is  shown  by  the  arrow  in  Figure  9.16  (the  same  as  Figure  6.17).  The 
experimental  2  is  shown  to  fall  in  the  cluster  of  model-based  simulation  results  of 
the  convection  class.  The  classifier  thus  declares  that  the  AUV-measured  WEarth  at 
the  250-m  depth  in  the  mixed  layer  is  convective.  This  result  is  consistent  with 
what  meteorological  and  hydrographic  conditions  indicate  in  Section  9.2,  and  with 
the  Lagrangian  float  observations  given  in  Section  9.3. 

9.5.2  At  20-m  depth  of  AUV  Mission  B980410T 

The  PSD  estimate  of  WEarth  on  the  first  and  second  legs  at  the  20-m  depth  is  shown 
in  Figure  9.17,  using  five-point  frequency-domain  smoothing  to  reduce  the  variance 
of  estimate.  We  note  a  spectral  peak  at  0.002  Hz,  and  a  lower  spectral  peak  at 
0.007  Hz  (side-lobe  possibility  is  ruled  out  considering  the  data  length  of  1400  s  and 
the  smoothing  window  of  50  s).  Compared  with  the  counterpart  PSD  of  the  250-m 
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PSD  templates  versus  PSD' of  AUV  experimental  data  at  250-m  depth,  aN  at  AUV  speed  1  m/s 


Figure  9.15:  PSD  estimate  of  WEarth  at  the  250-m  depth  of  Mission  B9804107  (non- 
smoothed)  along  with  PSD  templates  of  internal  wave  and  convection. 

depth  (for  the  first  and  second  legs)  as  shown  in  Figure  9.14,  we  notice  that  the  two 
PSDs  both  have  spectral  peaks  at  0.007  Hz.  In  the  time  domain,  Figure  9.18  shows 
the  cross-covariance  function  of  w  Earth  between  the  two  depths  (on  the  2nd  leg  at  each 
depth).  The  function  has  a  period  corresponding. to  0.007  Hz,  and  has  an  amplitude 
of  about  20%  or  larger  (note  that  the  tapering  at  tails  is  largely  due  to  sliding  of  time 
series  for  covariance  calculation).  This  time-domain  observation  shows  consistency 
with  the  peak  frequency  of  0.007  Hz  that  is  shared  by  the  two  depths. 

The  convection  PSD  template  used  by  the  AUV-based  classifier  is  developed  from 
the  MIT  Convection  Model  data.  The  model  uses  a  mixed  layer  of  350-m  depth  which 
is  set  to  the  same  value  as  that  measured  by  the  AUV  during  Mission  B9804107.  Due 
to  its  0.007-Hz  peak  frequency  which  is  close  to  the  model  template’s  0.005  Hz,  the 
250-m  depth  vertical  flow  velocity  data  is  determined  to  be  convective  by  the  classifier. 
We  thus  have  reason  to  believe  that  the  0.007-Hz  spectral  peak  is  associated  with 
ocean  convection. 

The  20-m  data,  however,  has  an  additional  higher  peak  at  0.002  Hz.  Due  to  this 
low-frequency  peak,  the  20-m  depth  vertical  flow  velocity  differs  from  the  model- 
based  convection  template,  thus  making  classification  difficult.  The  classification  test 
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Figure  9.16:  At  AUV  speed  1  m/s,  histograms  of  feature  z  for  internal  wave  and 
convection.  The  value  of  2  corresponding  to  the  250-m  data  of  Mission  B9804107  is 
marked  by  the  arrow. 

result  for  the  20-m  depth  wEarth  is  shown  in  Figure  9.19,  marked  by  the  arrow.  The 
arrow  still  lies  in  the  convection  cluster,  but  at  the  very  tail  towards  the  internal 
wave  cluster.  In  Section  9.6,  we  will  propose  the  cause  of  the  discrepancy  between 
shallow- water  data  and  the  ocean  convection  model. 


9.6  Discussions 

9.6.1  Effects  of  Ocean- Atmospheric  Coupling  on  Shallow- Water 
Measurement 

The  sea  water  loses  heat  to  the  air,  thus  convection  can  also  happen  in  the  atmosphere. 
Because  of  its  closeness  to  the  air-sea  boundary,  water  at  the  20-m  depth  is  much  more 
influenced  by  the  atmospheric  forcing  than  at  the  250-m  depth.  Vertical  flow  velocity 
at  the  20-m  depth  is  therefore  expected  to  bear  characteristics  of  the  atmospheric 
structure.  Ocean-atmosphere  coupling  is  not  included  in  the  convection  model.  It  is 
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PSD  of  the  1st  and  2nd  legs,  20-m  depth.  After  5-point  freq.  smoothing  with  Af=7  x  ICf4  Hz. 


Figure  9.17:  PSD  estimate  of  w^arth  at  the  20-m  depth  of  AUV  Mission  B9804107. 
The  l-cr  error  band  is  shown,  using  five-point  frequency-domain  smoothing. 

probably  this  inadequacy  that  has  caused  the  discrepancy  between  the  20-m  depth 
AUV  data  and  the  ocean  convection  model,  and  consequently  a  difficult  classification. 

The  atmospheric  spatial  scale  is  not  the  same  as  the  ocean’s,  and  is  probably 
larger  due  to  a  thicker  convective  layer.  Hence  it  is  quite  possible  that  the  low- 
frequency  (0.002-Hz)  spectral  peak4  (which  makes  classification  difficult)  of  the  20-m 
depth  data  is  due  to  the  atmospheric  forcing.  In  modeling  ocean  processes  in  shallow 
water,  it  is  thus  necessary  to  include  ocean-atmosphere  coupling  [86],  [87],  [88].  Good 
performance  of  AUV-based  classification  relies  on  well  modeled  templates. 

The  20-m  depth  data  also  brings  out  the  issue  of  adjusting  the  AUV  speed.  Ac¬ 
cording  to  the  illustration  in  Figure  6.4,  the  mingled  spectrum  peak  is  projected  from 
the  convection’s  spatial  periodicity.  The  peak  frequency  in  the  mingled  spectrum  in¬ 
creases  with  the  vehicle  speed.  At  the  20-m  depth,  convection’s  spatial  scale  is  larger 
than  that  in  deeper  water,  probably  due  to  the  atmospheric  forcing.  At  an  AUV 
speed  of  1  m/s,  this  larger  spatial  scale  is  projected  to  the  0.002-Hz  low-frequency 

4In  consideration  of  the  typical  surface  wave  period  of  10  s  [59],  it  is  not  very  likely  that  surface 
waves  are  the  cause  of  this  low-frequency  peak. 
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Cross-covariance  function  of  w£arth  between  250-m  and  20-m  depths  (2nd  leg  of  each  depth) 


Figure  9.18:  Cross-covariance  function  of  WEarth  between  250-m  and  20-m  depths  (on 
the  2nd  leg  at  each  depth) . 

peak  on  the  mingled  spectrum.  If  we  increase  the  vehicle  speed,  this  frequency  of 
convection  will  be  higher,  offering  a  better  separation  from  internal  wave’s  mingled 
spectrum  which  remains  at  the  base  band.  Therefore,  we  should  increase  the  AUV 
speed  when  the  spatial  scale  of  convection  is  larger. 

9.6.2  Unsteadiness  of  Surface  Cooling 

Let  us  compare  the  root-mean-square  (rms)  values  of  vertical  flow  velocity  given  by 
the  model  and  shown  by  the  Labrador  Sea  AUV  data  during  Mission  B9804107: 
wAUV= 0.020  m/s  <  wmodei= 0.043  m/s  at  the  250-m  depth 
*F4C/v=0-035  m/s  >  wmodei= 0.014  m/s  at  the  20-m  depth 

The  two  pairs  of  values  agree  within  a  factor  of  less  than  2.5.  What  attracts  our 
attention,  however,  is  the  relative  magnitude  of  wAuv  at  the  20-m  depth.  It  is  not 
only  larger  than  the  modeled  counterpart,  but  also  larger  than  wAuv  at  the  250-m 
depth.  One  possible  explanation  is  associated  with  the  ocean-atmospheric  coupling 
which  is  not  included  in  the  model.  The  coupling  may  enhance  convective  vertical 
velocity  in  shallow  water. 

Stronger  sea  surface  cooling  would  also  enhance  convection.  According  to  the 
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Class  1:  internal  wave.  AUV  speeds 1  m/s.  m^-1.15,  oz1=2.08 
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Figure  9.19:  At  AUV  speed  1  m/s,  histograms  of  feature  z  for  internal  wave  and 
convection.  The  value  of  z  corresponding  to  the  20-m  depth  data  of  Mission  B9804107 
is  marked  by  the  arrow. 


convection  model  theory  [53],  flow  velocity  amplitude  is  proportional  to  the  square 
root  of  surface  heat  flux.  AUV  runs  at  the  two  depths  were  both  conducted  shortly 
after  midnight  (local  time).  The  20-m  depth  run  was  about  one  hour  later  than 
the  250-m  depth  run.  Stronger  night  cooling  might  have  contributed  to  the  larger 
convective  vertical  velocity  at  the  20-m  depth.  In  the  model  setting  presented  in  Sub¬ 
section  5.1.3,  we  apply  a  steady  surface  heat  flux  of  300  W/m2  which  is  an  average 
value  around  the  AUV  mission  time.  For  more  accurate  description,  unsteadiness  of 
surface  cooling  should  be  incorporated  when  setting  the  model  parameters. 

Consideration  of  unsteady  atmospheric  forcing  may  also  help  us  understand  the 
heat  flux  calculation  results  based  on  the  AUV’s  measurements:  Earth-referenced 
vertical  flow  velocity  w Earth  and  potential  temperature  T.  The  average  heat  flux  in 
duration  r  can  be  calculated  by  [9] 


Q  = 


(9-1) 


where  po=1028  kg/m3  is  the  water  density  and  (7^=3900  J/(kg  °C)  is  the  specific  heat 
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capacity  of  water.  w'Earth (t)  and  T'(t)  are  deviations  from  their  respective  means  in 
duration  r: 


Earth  (^) 


i  r 

^Earthij)  ~  I  ^Earth(j)dt 

T  Jo 


(9.2) 


T{t)  =  T(t)  -  -  [  T(t)dt  (9.3) 

T  Jo 

Applying  the  above  computation  method  to  the  whole  duration  of  the  20-m  depth 
run  (four  legs  each  lasting  for  720  s),  the  average  heat  flux  is  130  W/m2.  Note  that  in 
the  thesis  definition,  positive  WEarth  points  upward,  so  heat  flux  upward  should  have 
a  “+”  sign.  Therefore,  the  heat  flux  calculation  result  using  the  20-m  depth  data 
has  a  sign  consistent  with  that  of  surface  cooling.  The  value  is  lower  than  300  W/m2 
estimated  from  meteorological  data,  but  within  a  factor  of  less  than  3. 

The  calculation  results  using  the  250-m  depth  data  are  complicated.  On  the  first 
and  second  legs  (720  s  each)  before  the  significant  temperature  dip  at  3200  s  (please 
see  the  fourth  panel  in  Figure  9.10),  the  average  heat  flux  is  140  W/m2.  The  sign  is 
consistent  with  that  of  surface  cooling,  and  the  value  is  close  to  the  calculated  average 
heat  flux  at  the  20-m  depth.  On  the  third  and  fourth  legs,  plus  the  temperature  dip 
which  connects  leg  2  and  leg  3,  the  calculated  average  heat  flux  is  -930  W/m2.  The 
sign  is  inconsistent  with  that  of  surface  cooling.  This  phenomenon  needs  further  study, 
but  we  can  propose  possible  explanations.  At  ~3  cm/s  convective  vertical  velocity, 
temperature  at  the  250-m  depth  is  the  result  of  more  than  two  hours  of  convection 
from  the  surface.  This  temperature  may  not  represent  the  present  conditions  at  the 
surface  if  the  atmospheric  forcing  is  unsteady.  If  the  surface  cooling  has  turned  off, 
there  could  be  rebound  of  convective  column  resulting  in  cold  water  going  up.  We 
should  note  that  simple  correlation  depends  on  a  steady-state  process. 

On  the  other  hand,  we  observe  that  the  time  series  of  WEarth  on  the  third  and  fourth 
legs  at  the  250-m  depth  no  longer  has  the  periodicity  shown  in  the  first  and  second 
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legs.  It  is  likely  that  due  to  unsteadiness  of  meteorological  conditions,  processes 
occurring  on  these  two  legs  differ  from  convection.  Besides  the  above  consideration  of 
unsteadiness,  a  closely  related  issue  is  the  length  of  AUV  surveys,  as  discussed  below. 

9.6.3  Consideration  of  AUV  Survey  Length 

Towards  the  end  of  Section  6.3,  we  relate  data  length  to  ueff  in  periodogram  estima¬ 
tion.  When  data  are  short,  we  have  to  tolerate  periodogram’s  variance  to  leave  enough 
uncorrelated  frequency  points  for  the  classifier’s  operation.  Another  issue  concerning 
data  length  is  ergodicity  [42],  [43],  [24],  [44],  In  the  time  domain,  an  ergodic  random 
process  is  defined  as  one  whose  first  and  second  moments  (ensemble  averages)  can 
be  replaced  by  the  corresponding  time  averages  as  the  data  length  approaches  infin¬ 
ity  [43].  Ergodicity  is  often  assumed  in  practical  data  analysis,  since  we  rarely  have 
an  ensemble  of  realizations  but  need  to  deal  with  a  single  realization  of  the  random 
process  [7].  A  longer  sample  time  series  would  better  cover  the  statistical  variations 
of  the  overall  random  process. 

The  concept  of  ergodicity  is  equally  applicable  in  the  space  domain.  A  longer 
AUV  survey  would  better  cover  the  statistical  variations  of  the  random  spatial  field. 
In  AUV  Mission  B9804107,  each  leg  had  a  length  of  only  about  800  m.  By  inspecting 
the  upper  panel  in  Figure  5.2,  we  notice  that  over  a  scale  of  up  to  several  hundred 
meters,  an  AUV  may  encounter  a  continuous  “row”  of  positive  or  negative  vertical  flow 
velocity.  This  localized  inhomogeneity  can  be  overcome  by  a  long  AUV  run  since  the 
orientation  of  such  rows  is  random.  For  an  AUV  leg  shorter  than  1  km  as  in  Mission 
B9804107,  ergodicity  of  the  measurement  would  be  challenged  if  the  AUV  happened 
to  fly  along  a  row.  As  shown  in  the  first  panel  in  Figure  9.11,  wEarth  has  a  negative 
mean  on  the  third  and  fourth  legs  at  the  250-m  depth  during  Mission  B9804107. 
While  keeping  in  mind  other  considerations  (measurement  bias  in  Section  8.4  and 
unsteadiness  of  surface  cooling  in  Subsection  9.6.2)  regarding  these  two  legs,  we  should 
also  consider  the  possibility  that  the  AUV  ran  along  downwelling  rows. 

The  shortness  of  AUV  runs  also  makes  the  heat  flux  calculation  difficult.  The 
mean  values  of  vertical  flow  velocity  and  temperature  are  removed  in  applying  Equa- 
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tion  (9.1).  These  mean  values  are  for  the  data  duration  r,  as  shown  in  Equation  (9.2) 
and  Equation  (9.3).  A  reliable  heat  flux  estimation  requires  long  observations,  i.e., 
long  AUV  runs. 

From  the  above  two  viewpoints,  AUV  survey  lines  should  be  made  longer.  This 
lengthening,  however,  is  not  without  limitation.  As  we  have  discussed  in  Subsec¬ 
tion  9.6.2,  the  assumption  of  stationarity  of  the  field  would  be  compromised  over  a 
longer  duration  due  to  possible  unsteadiness  of  the  atmospheric  forcing.  In  AUV 
survey  designs,  the  length  of  cruise  legs  should  be  determined  by  trading  off  the 
requirement  of  spatial  ergodicity  against  that  of  temporal  stationarity. 
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Chapter  10 


Conclusions  and  Future  Work 

10.1  Conclusions 

1.  The  “mingled  spectrum  principle”  concisely  relates  observations  from  a  moving 
platform  to  the  temporal-spatial  spectrum  of  the  process  under  survey. 

2.  The  principle  can  be  utilized  to  optimize  surveys  for  classifying  ocean  processes. 
It  lays  the  theoretical  basis  for  designing  an  AUV-based  classifier. 

3.  AUV-based  classification  is  demonstrated  for  distinguishing  ocean  convection 
from  internal  waves.  Simulation  results  show  that  at  a  higher  vehicle  speed,  the 
classification  performance  is  better  as  the  distinction  between  convection  and 
internal  wave  is  highlighted. 

4.  Tests  using  model  data  and  field  data  demonstrate  that  we  can  utilize  the  AUV’s 
controllable  motion  to  the  advantage  of  ocean  process  classification. 

10.2  Summary  of  Contributions 

1 .  We  established  the  “mingled  spectrum  principle” . 

2.  By  utilizing  the  mingled  spectrum  principle,  we  developed  a  parametric  tool 
for  designing  an  AUV-based  spectral  classifier.  Based  on  general  methods  of 
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feature  extraction,  we  presented  an  approach  for  classifying  Power  Spectrum 
Density  (PSD). 

3.  In  consideration  of  model  parameter  uncertainties,  we  devised  a  method  to  make 
the  classifier  robust  to  mismatch  between  models  and  real  data. 

4.  We  applied  the  AUV-based  classifier  to  distinguish  ocean  convection  from  in¬ 
ternal  wave. 

5.  We  installed  a  high-precision  acoustic  Doppler  sonar  in  an  AUV  to  measure  flow 
velocity.  We  calibrated  this  AUV-borne  measurement  system  in  a  high-precision 
tow  tank.  The  system  acquired  field  data  from  the  Labrador  Sea  in  1998.  We 
developed  the  data  processing  method  to  extract  the  Earth-referenced  vertical 
flow  velocity  from  AUV’s  raw  measurement,  and  the  error  analysis  approach. 

6.  Using  the  Labrador  Sea  flow  velocity  data  acquired  by  the  AUV,  the  classifi¬ 
cation  test  result  detects  convection’s  occurrence.  This  finding  is  supported  by 
more  traditional  oceanographic  analyses  and  observations. 

10.3  Future  Work 

10.3.1  Expand  Dimension  of  Classification  Quantities 

Classification  of  convection  versus  internal  wave  is  studied  as  an  example  application 
in  the  thesis.  Since  the  vertical  flow  velocity  is  a  key  signature  of  both  processes  [47], 
[16],  [2],  we  use  it  as  the  classification  quantity.  The  ocean  process  X ( t ,  r)  in  Figure  3.1 
and  Figure  4.1  thus  refers  to  the  vertical  flow  velocity.  Corresponding  to  this  single 
classification  quantity,  the  process  A  is  a  scalar  (as  a  function  of  time  and  space), 
i.e.,  of  one  dimension. 

To  fully  utilize  the  information  resources,  we  can  add  in  more  classification  quan¬ 
tities,  such  as  temperature.  The  addition  is  equivalent  to  expanding  the  dimension 
of  process  X.  For  example,  when  using  vertical  flow  velocity  and  temperature  as  the 
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two  quantities  for  classification,  X  will  become  a  two-dimensional  vector  with  each 
component  being  a  function  of  time  and  space. 

With  the  dimension  expansion,  the  mingled  spectrum  computation  should  be  cor¬ 
respondingly  extended.  Not  only  each  component’s  mingled  spectrum,  but  also  the 
cross-mingled-spectra  between  components,  will  be  useful  for  classification.  The  cross¬ 
spectrum  between  vertical  flow  velocity  and  temperature,  for  instance,  will  reflect  the 
correlation  between  the  two  quantities.  We  should  note,  however,  this  correlation  is 
based  on  the  AUV’s  “mingled”  measurements.  As  we  have  demonstrated,  the  vehicle 
speed  is  still  the  tuning  factor  we  should  optimize  for  good  classification. 

10.3.2  Extend  Two-Class  Method  to  Multi-Class  Method 

We  consider  a  two-class  problem  in  the  thesis.  In  distinguishing  convection  from 
internal  wave,  the  binary  hypothesis  formulation  is  plausible,  as  the  former  process 
occurs  in  a  vertically  mixed  water  column  while  the  latter  occurs  in  a  stably  stratified 
water  column.  To  treat  more  complicated  ocean  process  possibilities,  we  should 
extend  the  binary  classification  method  to  M-ary  [21]  (i.e.,  multi-class)  classification. 

The  class  separability  metric  (Equation  (2.24))  based  on  the  within-class  scatter 
matrix  Aw_y  (Equation  (2.22))  and  the  between-class  scatter  matrix  Ab_y  (Equa¬ 
tion  (2.23))  is  readily  applicable  to  L-class  problems  [3],  noting  the  definition  of  the 
overall  mean  M0  =  where  L  ^  2.  Correspondingly,  the  feature  extrac¬ 

tion  will  map  the  observation  vector  to  an  (L-l)-dimensional  feature  vector  Z.  For  a 
two-class  problem  where  L  =  2,  vector  Z  reduces  to  a  scalar  feature  z,  as  seen  in  the 
thesis. 

10.3.3  Incorporate  Real-Time  Algorithm  into  AUV  Software 

To  ensure  quick  response  from  the  AUV,  the  classification  algorithm  should  be  im¬ 
plemented  in  real  time.  In  the  thesis  research,  we  have  made  efforts  to  meet  this 
requirement: 
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1.  The  AUV-based  spectral  classifier  is  linear.  As  we  have  commented  at  the  end 
of  Section  2.3,  light  computational  load  is  a  major  reason  for  our  formulating  a 
linear  classifier. 

2.  The  classifier  has  demonstrated  its  performance  using  a  short  data  segment,  as 
presented  in  Section  6.5,  Section  7.4,  and  Section  9.5.  Capability  of  classification 
within  a  short  duration  is  important  for  real-time  or  quasi-real-time  operations, 
as  discussed  at  the  end  of  Section  6.3. 

To  make  the  classification  algorithm  operational  in  an  AUV,  we  need  to  eventually 
integrate  it  into  the  vehicle’s  software  architecture. 
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Appendix  A 


Relations  between  Various  Forms 
of  Temporal-Spatial  PSD 


By  definition  [89],  [25],  the  “complete”  autocorrelation  function  of  a  stationary  and 
homogeneous  field  X(t,x i,x2)  is 


Rxjid {r,pi,p2)  =  E[X(t,xi,x2)  ■  X(t  +  r,xi+  pi,x2  +  p2)}  (A.l) 


and  its  Fourier  transform,  i.e.,  the  “complete”  Power  Spectrum  Density  (PSD)  is 


/oo  roo  /*( 

/  / 

oo  J  — oo  J  — oo 


R*_3d(t,  p1,p2)e-j2^Tej2^eR^drdp1dp2 

(A.2) 


Now  let  us  integrate  Sxj3d{v^  vi,  ^2)  over  u2\ 
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/oo 

SxJ3d(Vi  uli  Vi)dV2 

-OO 

/oo  poo  poo  poo 

/  /  /  RxjsD(T,p1,P*)e-l2”rrei2’™e’2™dTdp1dp2dv2 

•oo  J  — oo  «/  — oo  J  —oo 

/OO  /‘OO  /‘OO 

/  /  RxjtD(T,PLP2)e-*'"1Tet2™>’'6(p2)dTdpIdp2  (A3) 

/OO  /'OO 

/  p1,0)e-*™eP™''dTdpl 

-oo  J  — oo 

On  the  other  hand,  for  two  locations  on  a  line  along  the  aq-axis,  the  “line”  auto¬ 
correlation  function  is 


Rx(t, Pi)  =  E[X(t,x i,x2)  •  AT(t  +  r,  aq  +  pi,x2)]  =  -R*-3£>(dPi,  0) 

(A.4) 


Then,  replacing  Rx_3D(t,  pi,  0)  with  #x(r,pi)  in  Equation  (A. 3),  we  get 


/OO 

SxJ3d(Vi  ^1)  u2)dU2 

-OO 

/OO  /‘OO 

/  Rx(r.Pi)e-j2^Tej2^dTdPl 
-OO  J  —oo 

=  Sxfavi) 


(A.5) 


where  the  last  step  is  by  the  definition  that  the  “line”  PSD  is  just  the  Fourier  trans¬ 
form  of  the  “line”  autocorrelation  function.  So  we  see  that  the  “line”  PSD  equals 
the  integration  of  the  “complete”  PSD  over  one  spatial  frequency.  Equation  (A.5) 
generally  holds  no  matter  whether  the  field  is  isotropic.  The  notation  “Sx{v,  v)n  used 
in  the  derivation  of  the  mingled  spectrum  formula  in  Subsection  3.1.1  means  the  same 
as  aSx{rhl'iT  in  Equation  (A.5). 

If  the  field  is  furthermore  isotropic,  its  spatial  spectrum  will  also  be  isotropic  [90]. 
Equation  (A.5)  is  then  illustrated  in  Figure  A.l.  For  an  isotropic  spectrum,  only  one 
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Figure  A.l:  Obtaining  the  “line”  PSD  by  integrating  the  “complete”  PSD  (contours 
shown  are  for  an  isotropic  field). 

spatial  frequency  suffices  for  description.  Thus  the  “complete”  PSD  can  be  written 
as  a  “radial”  PSD:  Sx_jrad.ia.iiji ^radial)-  Note  that  at  any  spatial  frequency  i 'radiai, 
the  “radial”  PSD  implicitly  represents  the  integration  on  a  circle  of  radius  2nvTadiai. 
Hence  the  “point-wise”  PSD  is 


S X -point— wiseiV j  ^ radial ) 


1 


2vr  I'radial 


Sx -radial  {V:  ^ radial ) 


(A.6) 


Based  on  isotropy,  we  know  that 


Sx.3d(V,  Vu  ^2 )  S x -point— wiseiV >  A  ^2) 


(A- 7) 


Then  the  three  variants  of  the  “complete”  PSD  of  an  isotropic  field  are  related  by 
the  following  equation: 
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Sx-3d{V,  u2)  ~  Sx  -point— wise 

1 


‘I'KyJvl  +  v\ 


(v, 

Sx  -radial  (jh  \J v\  +  v\) 


(A-8) 


Incorporating  Equation  (A. 8)  into  Equation  (A. 5),  we  obtain  Sx(v,^ i)  from 
S x jradiaiijli  Vui  +  vl)  bY  a  simple  integration. 

For  any  real  process  X(t,xi)  (on  a  line  along  the  aq-axis),  we  have 


Rx{r,pi)  =  Rx{-T,-pi) 


(A-9) 


since  E[X(t,Xi)  ■  X (t  +  r,  aq  4-  pi)]  =  E[X(t  +  r,  Xi  +  pi)  ■  X(t,  aq)]. 

By  the  property  of  two-dimensional  Fourier  transform  [65],  we  obtain  the  general 
symmetry  of  the  “line”  PSD: 


Sx(’1,‘'i)  =  Sx(-’l,  -"i) 


(A-10) 


When  the  field  is  furthermore  isotropic,  we  will  additionally  have 


Rx{r,Pi)  =  Rx{r,~pi ) 


(All) 


since  the  autocorrelation  does  not  depend  on  the  spatial  direction.  It  leads  to1 


Sx{r),vi )  =  Sx{t},-Ui) 


(A.12) 


1Equation  (A.12)  can  be  derived  in  another  way:  Sx(rj, -V\)  =  Sx-Zd(jU  ~vu  V2)dv2  = 
fT00Sx-3D(r],vi,i'2)dv2  =  Sx{r],v  1)  where  SX-3d(.V,  ~v\,  v2)  =  SX-3d{v,Vi,v2)  is  due  to 
isotropy  [2], 
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Combining  Equation  (A. 10)  and  Equation  (A. 12),  we  see  that  for  an  isotrpic  field, 
Sx{f 7,  v i)  is  symmetric  about  the  r]- axis  and  also  about  the  zq-axis. 
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Appendix  B 


Mechanical  Design  for  the 
Calibration  Experiment 


In  Section  8.3,  we  have  presented  the  calibration  experiment  for  the  AUV-borne  flow 
measurement  system  in  the  David  Taylor  Model  Basin.  In  this  appendix,  we  give 
details  of  the  mechanical  design.  The  mechanical  system  is  composed  of  three  parts: 
a  bracket,  a  wedge,  and  a  hull  platform,  as  shown  in  Figure  B.l. 

1.  The  bracket  is  for  fixing  the  mechanical  system  to  the  strut  of  the  tow  tank 
carriage.  Its  upper  rectangular  plate  attaches  to  the  base  of  the  strut  via  six 
bolts.  Its  circular  lower  plate,  as  shown  in  Figure  B.2,  attaches  to  the  wedge 
via  four  bolts.  The  multiple  hole  positions  are  designed  for  adjustability  of  the 
vehicle’s  yaw  angle. 

2.  The  wedge  is  for  adjusting  the  vehicle’s  pitch/roll  angle.  Four  wedges  of  differ¬ 
ent  angles  are  built.  Their  side  views  and  views  perpendicular  to  their  upper 
intersections  are  shown  in  Figure  B.3  and  Figure  B.4,  respectively.  The  wedge’s 
upper  elliptical  surface  attaches  to  the  bracket,  and  its  lower  circular  surface 
attaches  to  the  hull  platform. 

3.  The  hull  platform  is  for  holding  both  the  ADV  probe  and  the  vehicle’s  hull. 
It  is  attached  to  the  wedge  via  its  upper  circular  plate.  The  ADV  probe  is 
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held  by  the  45°  slanted  arm  as  shown  in  Figure  B.5.  The  lower  rectangular 
plate  attaches  to  the  flat  section  of  the  vehicle’s  inner  fairing,  and  thus  upholds 
the  whole  vehicle.  To  achieve  rigid  attachment,  the  lower  plate  matches  the 
vehicle’s  inner  fairing  with  the  largest  possible  area.  The  slanted  arm  ensures 
that  the  ADV  probe’s  stem  lies  in  the  vehicle’s  largest  vertical  cross-section. 
The  farther  the  ADV’s  sampling  volume  from  the  hull,  the  smaller  influence 
of  the  vehicle’s  hull  on  the  flow.  To  minimize  the  influence,  the  ADV  probe  is 
positioned  such  that  its  three  receivers  reach  the  brink  of  the  hull  but  do  not 
protrude. 

A  photo  showing  the  system  is  displayed  in  the  lower  panel  of  Figure  8.5.  The 
vehicle  is  mounted  on  the  carriage  strut  using  a  5°  wedge.  The  tips  of  the  ADV’s 
three  receivers  are  visible. 
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Figure  B.l:  The  overall  structure  of  the  mechanical  system  for  the  calibration  exper¬ 
iment  (“nipple”  is  renamed  as  “bracket”  in  the  thesis). 
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Figure  B.2:  Top  view  of  the  lower  plate  of  the  bracket. 
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Figure  B.3:  Side  views  of  the  four  wedges. 
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Figure  B.5:  Top  view  and  side  view  of  the  45°  slanted  arm  for  holding  the  ADV. 
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16.  Abstract  (Limit:  200  words) 

This  thesis  develops  and  demonstrates  methods  of  classifying  ocean  processes  using  an  underwater  moving  platform 
such  as  an  Autonomous  Underwater  Vehicle  (AUV).  The  "mingled  spectrum  principle"  is  established  which  concisely 
relates  observations  from  a  moving  platform  to  the  frequency-wavenumber  spectrum  of  the  ocean  process.  For  classifying 
different  processes,  an  AUV  is  not  only  able  to  jointly  utilize  the  time-space  information,  but  also  at  a  tunable  proportion 
by  adjusting  its  cruise  speed.  Based  on  the  mingled  spectrum  principle,  a  parametric  tool  for  designing  an  AUV-based 
spectral  classifier  is  developed. 

As  a  case  study,  AUV-based  classification  is  applied  to  distinguish  ocean  convection  from  internal  waves.  To  allow  for 
mismatch  between  modeled  templates  and  real  measurements,  the  AUV-based  classifier  is  designed  to  be  robust  to 
parameter  uncertainties.  By  simulation  tests  on  the  classifier,  it  is  demonstrated  that  at  a  higher  AUV  speed,  convection's 
distinct  spatial  feature  is  highlighted  to  the  advantage  of  classification. 

Experimental  data  are  used  to  test  the  AUV-based  classifier.  An  AUV-borne  flow  measurement  system  is  designed  and 
built  using,  an  Acoustic  Doppler  Velocimeter  (ADV).  In  February  1998,  the  AUV  acquired  field  data  of  flow  velocity  in  the 
Labrador  Sea  Convection  Experiment.  The  classification  test  result  detects  convection's  occurrence.  The  thesis  work 
provides  an  important  foundation  for  future  work  in  autonomous  detection  and  sampling  of  oceanographic  processes. 
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